ITS2_Caulosphere_Microbiome2017_FA19
R Packages
The following code uses several packages including vegan for community ecology, tidyr and plyr for organizing data tables and manipulating strings of text, and ggplot for plotting data.
#In this chunk I load all of the requisite packages.
library(vegan)
library(ggplot2)
library(ggpubr)
library(tidyr)
library(plyr)
library(data.table)
library(RColorBrewer)
library(ape)
library(car)
library(Hmisc)
library(plotly)## Warning: package 'plotly' was built under R version 3.6.2Data Import and Subsetting
Data Import
Below ITS data from drill shavings of phloem tissues collected from 47 Juglans nigra trees are imported into R. This includes a rarefied OTU table from mothur, the sample metadata which indicates the state from which the sample was collected and clone ID, and the taxonomic classifications for each OTU. Samples were rarefied to 3400 sequences in mothur prior to import, resulting in the loss of one sample from WA. Taxonomic assignments were made in mothur using the UNITE database.
# Below I am importing the rarefied OTU table for ITS2 from the caulosphere. This table was rarefied in Mothur by Geoff prior to importing into R to 3400 sequences per sample.
its2.ds.rareotu.wsingles<-read.table("Mothur_output/barkits2otutable.rarefied_jnigra17_fa19_gw.shared", header=TRUE, sep="\t",row.names = 2)
#Next I import the metadata for the particular habitat of interest. In this case I am working with drill shavings. I will use this metadata to subset the larger OTU file above.
its2.ds.met<-read.table("Mothur_output/barkits2metadata_jnigra17_su19_ajo.txt",header=TRUE,sep="\t",fill=TRUE,strip.white=TRUE)
#Next I import the taxonomy file.
its2.taxonomy.ds<-read.table("Mothur_output/barkits2taxonomy_jnigra17_fa19_gw.cons.taxonomy",header=TRUE,sep="\t",fill=TRUE,strip.white=TRUE,row.names=1)
#I then create a new object that contains the group names from the metadata file.
its2.ds.names<-its2.ds.met$GroupData Subsetting
Below I subset the OTU table to remove Mothur inserted metadata which indicates the percent identity threshold and the number of OTUs. I then remove singleton OTUs, OTUs that have only one 1 sequence following rarefaction.
#Below I remove the mothur associated metadata present in the first three columns to generate a matrix that consists of only OTU raw counts.
its2.ds.rareotu<-as.data.frame(its2.ds.rareotu.wsingles[,3:2696])
#I am now removing singleton OTUs (OTUs with only 1 representative sequence following rarefaction). This is to limit the effects of rare species on our distance matrices used during ordination.
its2.ds.rareotu.nosingletons<-as.data.frame(its2.ds.rareotu[,colSums(its2.ds.rareotu)>1])
#The following line is for saving the rarefied OTU table with singletons removed. This file is later reloaded to be used in downstream analyses. This is done to maintain consistency of results because rarefying generates slightly different results each time.
#write.table(its2.ds.rareotu.nosingletons, "~/Documents/microbiome2017/its2_bark_otutable_rarefied_nosingleton_fa19_ajo.shared",sep="\t", row.names=TRUE)Good’s Coverage
Relative Abundance Charts
To visualize differences in the relative abundance of different taxa I reorganize the OTU table and taxonomy strings to create relative abundance stacked bar charts. Below, stacked bar charts are made for phylum, class, and order levels.
Taxonomy Subsetting
To begin, the taxonomy file is subsetted to include only taxonomy information for OTUs present in the rarefied OTU table with singletons removed. The OTU table is first transposed so that OTU IDs are the row names and sample names are the column names. Then the taxonomic assignments are added to the OTU table. Taxonomic strings are then split by semi-colon into individual columns by Phylum, Class, Order, Family, Genus, and Species.
#Using the rarefied OTU tables I can generate relative abundance charts for the different taxonomic levels.To begin I first reformat the OTU table so that I can do some filtering and add taxonomy information.
its2.ds.rareotu.nosingles<-read.table("Mothur_output/barkits2otutable.rarefied.nosingleton_jnigra17_otutable_fa19_gw_ajo.shared", header=TRUE, sep="\t")
#I now transpose my OTU table so that I can begin to add in taxonomy information.
t.its2.ds.rareotu.nosingles<-as.data.frame(t(its2.ds.rareotu.nosingles))
#The below function is asking specifically for the OTU IDs. These are the row names. I will use these to subset the taxonomy file.
t.its2.ds.rareotu.nosingleslabs<-labels(t.its2.ds.rareotu.nosingles)
its2.ds.taxrare<-subset(its2.taxonomy.ds, rownames(its2.taxonomy.ds) %in% t.its2.ds.rareotu.nosingleslabs[[1]])
#Now I am subsetting the new taxonomy file to include only the taxonomy information stored in column 2.
its2.ds.taxrareinfo<-(its2.ds.taxrare[,2])
#Below I am using the separate function from tidyr to split the taxonomy strings into columns by semi-colon so that I can rename the OTUs to give them more meaning for downstream analyses.
library(tidyr)
#The UNITE database only provides classifications up to species level, thus, new columns are created up to species level.
its2.ds.taxonomylabs<-c("kingdom","phylum","class","order","family","genus","species")
#Below I create a new data frame with a column containing taxonomy information.
t.its2.ds.rareotutabtax<-data.frame(taxonomy=its2.ds.taxrareinfo,t.its2.ds.rareotu.nosingles)
#I then separate the taxonomy strings into new columns labeled with the taxonomy labels saved in the taxonomylabs object. Strings are being separated based on the presence of semi-colons.
t.its2.ds.rareotutabtaxsep<-separate(t.its2.ds.rareotutabtax,into=its2.ds.taxonomylabs,col=taxonomy,sep=";")Relative Abundance OTU Table
I then create a relative abundance OTU table. This table contains taxonomic information and can be searched to look at specific OTUs more easily.
#Using the rarefied OTU tables I can generate relative abundance charts for the different taxonomic levels.To begin I first reformat the OTU table so that I can do some filtering and add taxonomy information.
its2.ds.rareotu.nosingles.table<-read.table("Mothur_output/barkits2otutable.rarefied.nosingleton_jnigra17_otutable_fa19_gw_ajo.shared", header=TRUE, sep="\t")
rownames(its2.ds.rareotu.nosingles.table)<-make.names(row.names(its2.ds.rareotu.nosingletons))
its2.ds.rareotu.nosingles.table<-decostand(its2.ds.rareotu.nosingles.table,method="total")
#I now transpose my OTU table so that I can begin to add in taxonomy information.
t.its2.ds.rareotu.nosingles.table<-as.data.frame(t(its2.ds.rareotu.nosingles.table))
#The below function is asking specifically for the OTU IDs. These are the row names. I will use these to subset the taxonomy file.
t.its2.ds.rareotu.nosingleslabs.table<-labels(t.its2.ds.rareotu.nosingles.table)
its2.ds.taxrare.table<-subset(its2.taxonomy.ds, rownames(its2.taxonomy.ds) %in% t.its2.ds.rareotu.nosingleslabs.table[[1]])
#Now I am subsetting the new taxonomy file to include only the taxonomy information stored in column 2.
its2.ds.taxrareinfo<-(its2.ds.taxrare[,2])
#Below I am using the separate function from tidyr to split the taxonomy strings into columns by semi-colon so that I can rename the OTUs to give them more meaning for downstream analyses.
library(tidyr)
#The UNITE database only provides classifications up to species level, thus, new columns are created up to species level.
its2.ds.taxonomylabs<-c("kingdom","phylum","class","order","family","genus","species")
#Below I create a new data frame with a column containing taxonomy information.
t.its2.ds.rareotutabtax.table<-data.frame(taxonomy=its2.ds.taxrareinfo,t.its2.ds.rareotu.nosingles.table)
#I then separate the taxonomy strings into new columns labeled with the taxonomy labels saved in the taxonomylabs object. Strings are being separated based on the presence of semi-colons.
t.its2.ds.rareotutabtaxsep.table<-separate(t.its2.ds.rareotutabtax.table,into=its2.ds.taxonomylabs,col=taxonomy,sep=";")
#I then generate the html OTU table using the datatable function from DT package.
library(DT)
datatable(t.its2.ds.rareotutabtaxsep.table,fillContainer = T, height = "500px")Phylum Relative Abundance Chart
Below I generate a phylum relative abundance plot. I first extract the phylum information from the OTU table with taxonomy information. I then clean the taxa names to remove assignment confidence values and the p__ designator that precedes phylum assignment. OTU raw abundance is then summed by phylum and state, phyla that represent less than 1% of the total community are grouped into other, and relative abundance is calculated. Stacked bar charts are made using the melt function and ggplot.
#Below I create a data frame that consists of phylum assignment and the rarefied OTU table.
t.its2.ds.raretabphy<-data.frame(phylum=t.its2.ds.rareotutabtaxsep$phylum,t.its2.ds.rareotu.nosingles)
#I then remove the assignment confidence values contained in parentheses found in each taxonomic string using the sub function and the phylum "p__" designator at the start of each taxon.
t.its2.ds.raretabphy$phylum<-sub("\\(.*)","",x=t.its2.ds.raretabphy$phylum)
t.its2.ds.raretabphy$phylum<-sub("p__","",x=t.its2.ds.raretabphy$phylum)
#The following code is for the generation of phylum relative abundance stacked bar charts. I first sum each OTU by the phylum it belongs to.
library(plyr)
t.its2.ds.rareotutabphy<-as.data.frame(ddply(t.its2.ds.raretabphy, .(phylum),colwise(sum)))
rownames(t.its2.ds.rareotutabphy)<-make.names(t.its2.ds.rareotutabphy$phylum)
t.its2.ds.rareotutabphy<-t.its2.ds.rareotutabphy[,2:47]
#I then transpose the data frame so that sample names are now row names and column names are OTU names.
its2.ds.rareotutabphy<-as.data.frame(t(t.its2.ds.rareotutabphy))
its2.ds.rareotutabphy.state<-data.frame(state=its2.ds.met$State,its2.ds.rareotutabphy)
#I then sum each OTU by state and subset to exclude sample names. This is because the following functions only work on matrices containing numeric data.
its2.ds.rareotutabphy.statesum<-as.data.frame(ddply(its2.ds.rareotutabphy.state, .(state),colwise(sum)))
its2.ds.phylumcols<-its2.ds.rareotutabphy.statesum[,2:5]
#I am now creating an other category. Other includes those OTUs belonging to a phylum that comprises less than 1% of the total community. This new object is referred to as phyothers.
its2.ds.phyoth<-its2.ds.phylumcols[,colSums(its2.ds.phylumcols)/sum(its2.ds.phylumcols)<=0.01]
its2.ds.phyothers<-rowSums(its2.ds.phyoth)
#I then create an object that contains all of the phyla that comprise more than 1% of the total community.
its2.ds.phyreg<-its2.ds.phylumcols[,colSums(its2.ds.phylumcols)/sum(its2.ds.phylumcols)>0.01]
#I then create a new dataframe containing the state information, the other column, and the remaining phyla.
its2.ds.phytot2<-data.frame(state=its2.ds.rareotutabphy.statesum$state,its2.ds.phyreg,Other=its2.ds.phyothers)
#These values are then converted to relative abundance using the decostand function from vegan.
library(vegan)
its2.ds.phyrelabund<-decostand(its2.ds.phytot2[,2:4],method="total")
its2.ds.phyrelabund<-data.frame(state=its2.ds.rareotutabphy.statesum$state,its2.ds.phyrelabund)
#Below I use the melt function from the data.table package to reformat the data into stacked bar graph format.
library(data.table)
library(ggplot2)
library(ggpubr)
library(RColorBrewer)
its2.ds.phyrelabundmelt<-melt(its2.ds.phyrelabund, id.vars="state", variable.name="Phylum")
its2.ds.colors.n.phylum <- length(unique(its2.ds.phyrelabundmelt[,'Phylum']))
#Below I generate the relative abundance bar charts in ggplot
its2.ds.phyrel<-ggplot(its2.ds.phyrelabundmelt, aes(x=state, y=value, fill=Phylum))+
geom_bar(stat="identity", show.legend=TRUE, color="black")+
scale_fill_manual(values=colorRampPalette(brewer.pal(9, 'Set1'))(its2.ds.colors.n.phylum)) +
xlab("State") +
ylab("Relative Abundance") +
theme(panel.border = element_blank(),panel.background=element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.text.x = element_text(size = 14), axis.ticks.x = element_blank(), axis.text.y = element_text(size = 14), axis.line.y.left = element_line(), axis.title.y = element_text(size = 14))+
scale_y_continuous(breaks=c(0,0.25,0.50,0.75,1),limits=c(0,1.05))+
geom_text(data=NULL,aes(x=0.51, y=1.05,label="Caulosphere: ITS2"),hjust=0,colour="black")
ggplotly(its2.ds.phyrel, tooltip = c("Phylum"))Class Relative Abundance Chart
Below I generate a class relative abundance plot. I first extract the class information from the OTU table that contains taxonomy information. I then clean the taxa names to remove assignment confidence values, the p__ designator that precedes phylum assignment, and c__ designator that precedes class assignment. OTU raw abundance is then summed by class and state, classes that represent less than 1% of the total community are grouped into other, and relative abundance is calculated. Stacked bar charts are made using the melt function and ggplot.
#The following chunk creates a relative abundance bar chart at the class level.
library(plyr)
t.its2.ds.raretabcls<-data.frame(class=t.its2.ds.rareotutabtaxsep$class,t.its2.ds.rareotu.nosingles)
#What I am doing below is cleaning up the class names. These include assignment confidence levels, p__ and c__
t.its2.ds.raretabcls$class<-sub("\\(.*)","",x=t.its2.ds.raretabcls$class)
t.its2.ds.raretabcls$class<-sub("c__","",x=t.its2.ds.raretabcls$class)
t.its2.ds.raretabcls$class<-sub("p__","",x=t.its2.ds.raretabcls$class)
#I then sum the OTUs by class and make the row names the class name.
t.its2.ds.rareotutabcls<-as.data.frame(ddply(t.its2.ds.raretabcls, .(class),colwise(sum)))
rownames(t.its2.ds.rareotutabcls)<-make.names(t.its2.ds.rareotutabcls$class)
t.its2.ds.rareotutabcls<-t.its2.ds.rareotutabcls[,2:47]
#I then transpose the data frame and sum each class by state.
its2.ds.rareotutabcls<-as.data.frame(t(t.its2.ds.rareotutabcls))
its2.ds.rareotutabcls.state<-data.frame(state=its2.ds.met$State,its2.ds.rareotutabcls)
its2.ds.rareotutabcls.statesum<-as.data.frame(ddply(its2.ds.rareotutabcls.state, .(state),colwise(sum)))
#At the class level, UNITE may classify things as phylum_unclassified. These don't provide a lot of information and thus are grouped into a new group called unknown.
its2.ds.clsunknown<-its2.ds.rareotutabcls.statesum[,grep("_unclassified",names(its2.ds.rareotutabcls.statesum))]
its2.ds.clsunknownsum<-rowSums(its2.ds.clsunknown)
#I then use the same grep function as above to pull out OTUs classified to the classlevel. By setting invert=TRUE I select items lacking the strings in the grep command.
its2.ds.clstot<-its2.ds.rareotutabcls.statesum[,grep("_unclassified",names(its2.ds.rareotutabcls.statesum), invert=TRUE)]
#I then subset the data so that only numerical data is present and creating a new class of objects, other, that includes OTUs from classes that represent 1% or less of the total community.
its2.ds.clscols<-its2.ds.clstot[,2:23]
its2.ds.clsoth<-its2.ds.clscols[,colSums(its2.ds.clscols)/sum(its2.ds.clscols)<=0.01]
its2.ds.clsothers<-rowSums(its2.ds.clsoth)
its2.ds.clsreg<-its2.ds.clscols[,colSums(its2.ds.clscols)/sum(its2.ds.clscols)>0.01]
#I then create a new data frame that contains all of the class data and determine relative abundance using the decostand function from vegan.
its2.ds.clstot2<-data.frame(state=its2.ds.clstot$state,its2.ds.clsreg,Other=its2.ds.clsothers,Unclassified=its2.ds.clsunknownsum)
library(vegan)
its2.ds.clsrelabund<-decostand(its2.ds.clstot2[,2:13],method="total")
its2.ds.clsrelabund<-data.frame(state=its2.ds.rareotutabcls.statesum$state,its2.ds.clsrelabund)
#Below I use the melt function from the data.table package to reformat the data into stacked bar graph format.
library(data.table)
library(ggplot2)
library(ggpubr)
library(RColorBrewer)
its2.ds.clsrelabundmelt<-melt(its2.ds.clsrelabund, id.vars="state", variable.name="class")
its2.ds.cls.colorCount<- length(unique(its2.ds.clsrelabundmelt[,'class']))
getPalette=colorRampPalette(brewer.pal(9,"Set1"))
#Below I generate the relative abundance bar charts in ggplot
its2.ds.clsrel<-ggplot(its2.ds.clsrelabundmelt, aes(x=state, y=value, fill=class))+
geom_bar(stat="identity",show.legend=TRUE,color="black")+
scale_fill_manual(values=getPalette(its2.ds.cls.colorCount))+
xlab("State") +
ylab("Relative Abundance") +
theme(panel.border = element_blank(),panel.background=element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.text.x = element_text(size = 14), axis.ticks.x = element_blank(), axis.text.y = element_text(size = 14), axis.line.y.left = element_line(), axis.title.y = element_text(size = 14))+
scale_y_continuous(breaks=c(0,0.25,0.50,0.75,1),limits=c(0,1.05))+
geom_text(data=NULL,aes(x=0.50, y=1.05,label="Caulosphere: ITS2"),hjust=0,colour="black")
ggplotly(its2.ds.clsrel, tooltip=c("class"))Order Relative Abundance Chart
Below I generate a order relative abundance plot. I first extract the order information from the OTU table that contains taxonomy information. I then clean the taxa names to remove assignment confidence values, the p__ designator that precedes phylum assignment, c__ designator that precedes class assignment, and the o__ designator that precedes order assignment. OTU raw abundance is then summed by order and state, orders that represent less than 1% of the total community are grouped into other, and relative abundance is calculated. Stacked bar charts are made using the melt function and ggplot.
#The following chunk creates a relative abundance bar chart at the order level.
#What I am doing below is cleaning up the order names. These include assignment confidence levels, p__, c__, and o__
library(plyr)
t.its2.ds.raretabord<-data.frame(order=t.its2.ds.rareotutabtaxsep$order,t.its2.ds.rareotu.nosingles)
t.its2.ds.raretabord$order<-sub("\\(.*)","",x=t.its2.ds.raretabord$order)
t.its2.ds.raretabord$order<-sub("o__","",x=t.its2.ds.raretabord$order)
t.its2.ds.raretabord$order<-sub("c__","",x=t.its2.ds.raretabord$order)
t.its2.ds.raretabord$order<-sub("p__","",x=t.its2.ds.raretabord$order)
#I then sum the OTUs by class and make the row names the order name.
t.its2.ds.rareotutabord<-as.data.frame(ddply(t.its2.ds.raretabord, .(order),colwise(sum)))
rownames(t.its2.ds.rareotutabord)<-make.names(t.its2.ds.rareotutabord$order)
t.its2.ds.rareotutabord<-t.its2.ds.rareotutabord[,2:47]
#I then transpose the data frame and sum each order by state.
its2.ds.rareotutabord<-as.data.frame(t(t.its2.ds.rareotutabord))
its2.ds.rareotutabord.state<-data.frame(state=its2.ds.met$State,its2.ds.rareotutabord)
its2.ds.rareotutabord.statesum<-as.data.frame(ddply(its2.ds.rareotutabord.state, .(state),colwise(sum)))
#At the order level, UNITE may classify things as phylum_unclassified or class_unclassified. These don't provide a lot of information and thus are grouped into a new group called unknown.
its2.ds.ordunknown<-its2.ds.rareotutabord.statesum[,grep("_unclassified",names(its2.ds.rareotutabord.statesum))]
#I then use the same grep function as above to pull out OTUs classified to the order level. By setting invert=TRUE I select items lacking the strings in the grep command.
its2.ds.ordunknownsum<-rowSums(its2.ds.ordunknown)
its2.ds.ordtot<-its2.ds.rareotutabord.statesum[,grep("_unclassified",names(its2.ds.rareotutabord.statesum), invert=TRUE)]
#I then subset the data so that only numerical data is present and create a new class of objects, other, that includes OTUs from orders that represent 1% or less of the total community.
its2.ds.ordcols<-its2.ds.ordtot[,2:73]
its2.ds.ordoth<-its2.ds.ordcols[,colSums(its2.ds.ordcols)/sum(its2.ds.ordcols)<=0.01]
its2.ds.ordothers<-rowSums(its2.ds.ordoth)
its2.ds.ordreg<-its2.ds.ordcols[,colSums(its2.ds.ordcols)/sum(its2.ds.ordcols)>0.01]
#I then create a new data frame that contains all of the order data and determine relative abundance using the decostand function from vegan.
library(vegan)
its2.ds.ordtot2<-data.frame(state=its2.ds.ordtot$state,its2.ds.ordreg,Other=its2.ds.ordothers, Unclassified=its2.ds.ordunknownsum)
its2.ds.ordrelabund<-decostand(its2.ds.ordtot2[,2:14],method="total")
its2.ds.ordrelabund<-data.frame(state=its2.ds.rareotutabord.statesum$state,its2.ds.ordrelabund)
#Below I use the melt function from the data.table package to reformat the data into stacked bar graph format.
library(data.table)
library(ggplot2)
library(ggpubr)
library(RColorBrewer)
its2.ds.ordrelabundmelt<-melt(its2.ds.ordrelabund, id.vars="state", variable.name="order")
its2.ds.ord.colorCount<- length(unique(its2.ds.ordrelabundmelt[,'order']))
getPalette=colorRampPalette(brewer.pal(9,"Set1"))
#Below I generate the relative abundance bar charts in ggplot
its2.ds.ordrel<-ggplot(its2.ds.ordrelabundmelt, aes(x=state, y=value, fill=order))+
geom_bar(stat="identity",show.legend=TRUE,color="black")+
scale_fill_manual(values=getPalette(its2.ds.ord.colorCount))+
xlab("State") +
ylab("Relative Abundance") +
theme(panel.border = element_blank(),panel.background=element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.text.x = element_text(size = 14), axis.ticks.x = element_blank(), axis.text.y = element_text(size = 14), axis.line.y.left = element_line(), axis.title.y = element_text(size = 14))+
scale_y_continuous(breaks=c(0,0.25,0.50,0.75,1),limits=c(0,1.05))+
geom_text(data=NULL,aes(x=0.5, y=1.05,label="Caulosphere: ITS2"),hjust=0,colour="black")
ggplotly(its2.ds.ordrel,tooltip = c("order"))OTU Renaming
To give OTUs more meaningful names, taxonomy strings can be subset to change the OTU ID to unique taxonomic ID. Similar to the relative abundance charts, the taxonomy file must first me added to a transposed OTU table, taxonomy strings need to be split into different classification levels, and the lowest level of classification can be used to replace the original OTU ID. For instance, if an OTU received Ascomycota as its highest level of classification, the OTU ID would be changed to Ascomycota1.
Taxonomy Filtering
Taxonomy is first filtered to include only those OTUs present in the rarefied OTU table with no singletons. Then the OTU table is transposed so that OTU ID is the row name and sample ID is the column name. Taxonomy strings are then added to the OTU table and strings are split by semi-colon.
#Below I am reformatting the OTU table so that I can begin to filter and merge the taxonomy information.
#I first transpose the rarefied OTU table with no singletons.
t.its2.ds.rareotu.nosingles<-as.data.frame(t(its2.ds.rareotu.nosingles))
#The below function is asking specifically for the OTU IDs. These are the row names. I will use these to subset the taxonomy file.
t.its2.ds.rareotu.nosingleslabs<-labels(t.its2.ds.rareotu.nosingles)
its2.ds.taxrare<-subset(its2.taxonomy.ds, rownames(its2.taxonomy.ds) %in% t.its2.ds.rareotu.nosingleslabs[[1]])
its2.ds.taxrareinfo<-(its2.ds.taxrare[,2])
#Below I am using the separate function from tidyr to split the taxonomy strings into columns by semi-colon so that I can rename the OTUs to give them more meaning for downstream analyses.
library(tidyr)
its2.ds.taxonomylabs<-c("kingdom","phylum","class","order","family","genus","species")
t.its2.ds.rareotutabtax<-data.frame(taxonomy=its2.ds.taxrareinfo,t.its2.ds.rareotu.nosingles)
t.its2.ds.rareotutabtaxsep<-separate(t.its2.ds.rareotutabtax,into=its2.ds.taxonomylabs,col=taxonomy,sep=";")Renaming of OTUs
OTUs are then renamed. Prior to official renaming, taxonomy strings are cleaned so that assignment confidence values and taxon leaders (e.g. p__, c__, etc.) are removed.
#Here I am renaming the OTUs using the make.names function. This will make the OTU name something more informative. I am using species to do this.
t.its2.ds.rareotutabtaxsep$species<-sub("\\(.*)","",x=t.its2.ds.rareotutabtaxsep$species)
t.its2.ds.rareotutabtaxsep$species<-sub("(_unclassified)","",x=t.its2.ds.rareotutabtaxsep$species)
t.its2.ds.rareotutabtaxsep$species<-sub("(p__)","",x=t.its2.ds.rareotutabtaxsep$species)
t.its2.ds.rareotutabtaxsep$species<-sub("(c__)","",x=t.its2.ds.rareotutabtaxsep$species)
t.its2.ds.rareotutabtaxsep$species<-sub("(o__)","",x=t.its2.ds.rareotutabtaxsep$species)
t.its2.ds.rareotutabtaxsep$species<-sub("(f__)","",x=t.its2.ds.rareotutabtaxsep$species)
t.its2.ds.rareotutabtaxsep$species<-sub("(g__)","",x=t.its2.ds.rareotutabtaxsep$species)
t.its2.ds.rareotutabtaxsep$species<-sub("(s__)","",x=t.its2.ds.rareotutabtaxsep$species)
rownames(t.its2.ds.rareotutabtaxsep)<-make.names(paste(rownames(t.its2.ds.rareotutabtaxsep),t.its2.ds.rareotutabtaxsep$species,sep="_"))
t.its2.ds.rareotutabtax.rename<-t.its2.ds.rareotutabtaxsep[,8:53]
its2.ds.rareotutabtax.rename<-as.data.frame(t(t.its2.ds.rareotutabtax.rename))Principal Coordinate Analysis (PCoA)
To visually assess differences in the community composition of drill shaving fungal communities by clone and state, I perform a principal coordinate analysis (PCoA).
Relative Abundance Calculation
Prior to conducting a PCoA I first convert my OTU table from raw OTU abundances to relative abundance.
Principal Coordinate Analysis
I then perform the PCoA using the pcoa function from the ape package. Prior to performing the PCoA, a distance matrix needs to be calculated for the OTU table. In this case, I use the bray-curtis method.
#Below I calculate the distance matrix using the bray-curtis method and then perform the principal coordinate analysis on the relativized OTU table.
its2.ds.bacdist<-vegdist(its2.ds.relabund, method="bray")
library(ape)
its2.ds.bac.pcoa<-pcoa(its2.ds.bacdist)
#biplot(its2.ds.bac.pcoa,plot.axes = c(1,2))Plotting of Principal Coordinate Analysis
Then to plot the PCoA, I use the ggplot function. To use ggplot I need to first pull out the site scores for each sample and generate ellipses.
#We need to first subset our metadata table to include only those samples that passed rarefaction.
its2.ds.met.rare<-subset(its2.ds.met,its2.ds.met$Group%in%labels(its2.ds.rareotu.wsingles)[[1]])
#Then we extract PCoA site scores (these will be the x and y coordinates for each sample)
its2.ds.bac.pcoavec<-as.data.frame(its2.ds.bac.pcoa$vectors)
its2.ds.bac.pcoasitescores<-data.frame(PC1=its2.ds.bac.pcoavec$Axis.1, PC2=its2.ds.bac.pcoavec$Axis.2)
#Create a new dataframe that includes the site scores from above with metadata from your study.
its2.ds.bac.pcoagraph<-data.frame(its2.ds.bac.pcoasitescores,PC1=its2.ds.bac.pcoasitescores$PC1, PC2=its2.ds.bac.pcoasitescores$PC2, State=its2.ds.met.rare$State, Clone=its2.ds.met.rare$Clone,group=its2.ds.met.rare$Group)
#This is where you make confidence ellipses. I don't know what all the code means. Just know where to plug in my objects.
its2.ds.bac.pcoaellipse<-ordiellipse(its2.ds.bac.pcoasitescores,its2.ds.bac.pcoagraph$State, display="sites", kind="sd", draw="none")
df_ell.bark <- data.frame()
for(g in levels(its2.ds.bac.pcoagraph$State)){
df_ell.bark <- rbind(df_ell.bark, cbind(as.data.frame(with(its2.ds.bac.pcoagraph[its2.ds.bac.pcoagraph$State==g,], vegan:::veganCovEllipse(its2.ds.bac.pcoaellipse[[g]]$cov,its2.ds.bac.pcoaellipse[[g]]$center,its2.ds.bac.pcoaellipse[[g]]$scale))) ,State=g))}
#Plot PCoA in ggplot
library(ggplot2)
its2.ds.pcoa<-ggplot(its2.ds.bac.pcoagraph, aes(PC1,PC2, colour=State))+
#geom_text(aes(label=group))+
geom_path(data=df_ell.bark, aes(x=PC1, y=PC2, colour=State),size=0.5, linetype=1)+
geom_point(aes(shape=Clone), size=3.5)+
xlab("PC1 (34.7%)")+
ylab("PC2 (16.0%)")+
theme(axis.title.x=element_text(size=14, face="bold"))+
theme(axis.title.y=element_text(size=14, face="bold"))+
theme(axis.text.x=element_text(size=12, face="bold"))+
theme(axis.text.y=element_text(size=12, face="bold"))+
scale_color_manual(values=c("#006600","#3399FF","#FF9900"))+
scale_shape_manual(values=c(16,10,8,7,0,2))+
scale_x_continuous(breaks=c(-0.60,-0.30,0,0.30,0.60),limits=c(-0.60,0.60))+
scale_y_continuous(breaks=c(-0.60,-0.30,0,0.30,0.60),limits=c(-.60,.60))+
theme(panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
panel.border=element_rect(colour="black", size=1, fill=NA))+
theme(panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
panel.background = element_blank(),
panel.border=element_rect(color="black", size=1, fill=NA))+
geom_text(data=NULL,aes(x=-0.60,y=0.60,label="Caulosphere: ITS2"),hjust=0,colour="black")
its2.ds.pcoa
PERMANOVA
To formally test how well state and clone explain differences in drill shaving fungal community composition, I perform a PERMANOVA using the adonis function from vegan.
#Below I perform permanovas by state and clone using the adonis function from vegan.
library(vegan)
its2.ds.bac.perm.state.inter<-adonis(its2.ds.relabund~its2.ds.met.rare$State*its2.ds.met.rare$Clone, method="bray",permutations=10000)
its2.ds.bac.perm.state.inter##
## Call:
## adonis(formula = its2.ds.relabund ~ its2.ds.met.rare$State * its2.ds.met.rare$Clone, permutations = 10000, method = "bray")
##
## Permutation: free
## Number of permutations: 10000
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model
## its2.ds.met.rare$State 2 7.5527 3.7764 21.1822
## its2.ds.met.rare$Clone 4 1.4147 0.3537 1.9838
## its2.ds.met.rare$State:its2.ds.met.rare$Clone 8 2.0791 0.2599 1.4578
## Residuals 31 5.5267 0.1783
## Total 45 16.5732
## R2 Pr(>F)
## its2.ds.met.rare$State 0.45572 9.999e-05 ***
## its2.ds.met.rare$Clone 0.08536 0.006199 **
## its2.ds.met.rare$State:its2.ds.met.rare$Clone 0.12545 0.036896 *
## Residuals 0.33347
## Total 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Richness and Diversity
Then to assess the alpha diversity of our samples, I calculate the observed richness and shannon diversity index using the vegan package. This is done using the rarefied OTU table that contains singletons. Singletons are included to keep sampling depth standard between samples.
Richness
Observed OTU richness is first calculated for each sample. Then due to the unbalanced nature of our study design, we use a type 3 ANOVA.
#In this chunk we calculate richness and create a data frame including richness with metadata
its2.ds.met.rare<-read.table("Mothur_output/barkits2metadata_jnigra17_su19_ajo.txt",header=TRUE,sep="\t",fill=TRUE,strip.white=TRUE)
#I first import my rarefied OTU table that contains singletons.
its2.ds.rareotu.wsingles<-read.table("Mothur_output/barkits2otutable.rarefied_jnigra17_fa19_gw.shared", header=TRUE, sep="\t")
#I then remove mothur associated metadata present in the first 3 columns.
its2.ds.rareotu<-as.data.frame(its2.ds.rareotu.wsingles[,4:2697])
#Below I calculate observed species richness using the specnumber function from vegan. I then create a new data frame to include sample metadata.
library(vegan)
its2.ds.bac.rich<-specnumber(its2.ds.rareotu)
its2.ds.richnesstabmet<-data.frame(State=its2.ds.met.rare$State, clone=its2.ds.met.rare$Clone, sample=its2.ds.met.rare$Group,Richness=its2.ds.bac.rich)
#I first test for a significant interaction between state and clone. Due to the unbalanced design we use a type 3 ANOVA.
library(car)
its2.ds.anova.typ3.rich<-aov((Richness)~State*clone, data=its2.ds.richnesstabmet)
its2.ds.typ3aov.rich<-Anova(its2.ds.anova.typ3.rich,type="III")
its2.ds.typ3aov.rich## Anova Table (Type III tests)
##
## Response: (Richness)
## Sum Sq Df F value Pr(>F)
## (Intercept) 73633 1 127.7164 1.597e-12 ***
## State 8946 2 7.7580 0.001855 **
## clone 4708 4 2.0414 0.112875
## State:clone 9259 8 2.0075 0.078727 .
## Residuals 17873 31
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1| Sum Sq | Df | F value | Pr(>F) | |
|---|---|---|---|---|
| (Intercept) | 73633.333 | 1 | 127.716438 | 0.0000000 |
| State | 8945.542 | 2 | 7.757986 | 0.0018547 |
| clone | 4707.714 | 4 | 2.041373 | 0.1128755 |
| State:clone | 9259.161 | 8 | 2.007493 | 0.0787269 |
| Residuals | 17872.667 | 31 | NA | NA |
#The next line generates plots to evaluate regression assumptions.
par(mfrow=c(2,2))
plot(its2.ds.anova.typ3.rich)
#There was no significant interaction between state and clone and thus we drop the interaction term.
its2.ds.anova.add.rich<-aov((Richness)~State+clone, data=its2.ds.richnesstabmet)
its2.ds.add.aov.rich<-Anova(its2.ds.anova.add.rich,type="III")
its2.ds.add.aov.rich## Anova Table (Type III tests)
##
## Response: (Richness)
## Sum Sq Df F value Pr(>F)
## (Intercept) 151788 1 218.1846 < 2.2e-16 ***
## State 81162 2 58.3322 1.897e-12 ***
## clone 2526 4 0.9079 0.4689
## Residuals 27132 39
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Perform Tukey test
its2.ds.add.aov.rich.tukey<-TukeyHSD(its2.ds.anova.add.rich)
its2.ds.add.aov.rich.tukey.table<-rbind(its2.ds.add.aov.rich.tukey$State,its2.ds.add.aov.rich.tukey$clone)
its2.ds.add.aov.rich.tukey## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = (Richness) ~ State + clone, data = its2.ds.richnesstabmet)
##
## $State
## diff lwr upr p adj
## TN-IN 56.63492 33.73605 79.53379 0.0000014
## WA-IN -43.64286 -67.93078 -19.35494 0.0002522
## WA-TN -100.27778 -123.17665 -77.37891 0.0000000
##
## $clone
## diff lwr upr p adj
## 132-130 9.649471 -26.99877 46.29771 0.9423236
## 272-130 22.093915 -14.55432 58.74215 0.4316962
## 55-130 4.205357 -33.50537 41.91608 0.9976622
## WT-130 5.378968 -29.04606 39.80399 0.9914252
## 272-132 12.444444 -23.10957 47.99846 0.8533570
## 55-132 -5.444114 -42.09235 31.20412 0.9929267
## WT-132 -4.270503 -37.52824 28.98723 0.9959590
## 55-272 -17.888558 -54.53680 18.75968 0.6339629
## WT-272 -16.714947 -49.97268 16.54279 0.6081015
## WT-55 1.173611 -33.25141 35.59864 0.9999786| diff | lwr | upr | p adj | |
|---|---|---|---|---|
| TN-IN | 56.634921 | 33.73605 | 79.53379 | 0.0000014 |
| WA-IN | -43.642857 | -67.93078 | -19.35494 | 0.0002522 |
| WA-TN | -100.277778 | -123.17665 | -77.37891 | 0.0000000 |
| 132-130 | 9.649471 | -26.99877 | 46.29771 | 0.9423236 |
| 272-130 | 22.093915 | -14.55432 | 58.74215 | 0.4316962 |
| 55-130 | 4.205357 | -33.50537 | 41.91608 | 0.9976622 |
| WT-130 | 5.378968 | -29.04606 | 39.80399 | 0.9914252 |
| 272-132 | 12.444444 | -23.10957 | 47.99846 | 0.8533570 |
| 55-132 | -5.444114 | -42.09235 | 31.20412 | 0.9929267 |
| WT-132 | -4.270503 | -37.52824 | 28.98723 | 0.9959590 |
| 55-272 | -17.888558 | -54.53680 | 18.75968 | 0.6339629 |
| WT-272 | -16.714947 | -49.97268 | 16.54279 | 0.6081015 |
| WT-55 | 1.173611 | -33.25141 | 35.59864 | 0.9999786 |
#The below two lines generate plots to evaluate regression assumptions.
par(mfrow=c(2,2))
plot(its2.ds.anova.add.rich)
#There was no significant interaction between state and clone and thus we drop the interaction term.
its2.ds.anova.sing.rich<-aov((Richness)~State, data=its2.ds.richnesstabmet)
its2.ds.sing.aov.rich<-Anova(its2.ds.anova.sing.rich,type="III")
its2.ds.sing.aov.rich## Anova Table (Type III tests)
##
## Response: (Richness)
## Sum Sq Df F value Pr(>F)
## (Intercept) 391114 1 567.055 < 2.2e-16 ***
## State 80775 2 58.555 5.303e-13 ***
## Residuals 29658 43
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1| Sum Sq | Df | F value | Pr(>F) | |
|---|---|---|---|---|
| (Intercept) | 391114.29 | 1 | 567.0554 | 0 |
| State | 80774.65 | 2 | 58.5554 | 0 |
| Residuals | 29658.33 | 43 | NA | NA |
#The below two lines generate plots to evaluate regression assumptions.
par(mfrow=c(2,2))
plot(its2.ds.anova.sing.rich)
#Perform Tukey test
its2.ds.anova.sing.rich.tukey<-TukeyHSD(its2.ds.anova.sing.rich)
its2.ds.anova.sing.rich.tukey.table<-rbind(its2.ds.anova.sing.rich.tukey$State)
its2.ds.anova.sing.rich.tukey.table## diff lwr upr p adj
## TN-IN 56.63492 33.91733 79.35251 9.127502e-07
## WA-IN -43.64286 -67.73850 -19.54721 2.061864e-04
## WA-TN -100.27778 -122.99537 -77.56019 1.424305e-12| diff | lwr | upr | p adj | |
|---|---|---|---|---|
| TN-IN | 56.63492 | 33.91733 | 79.35251 | 0.0000009 |
| WA-IN | -43.64286 | -67.73850 | -19.54721 | 0.0002062 |
| WA-TN | -100.27778 | -122.99537 | -77.56019 | 0.0000000 |
#I then create a box plot for the richness data by state using ggplot.
library(Hmisc)
library(ggpubr)
its2.ds.richness.summarized<-summarize(its2.ds.richnesstabmet$Richness,by=its2.ds.richnesstabmet$State, FUN=max)
its2.ds.richness.summarized<-data.frame(State=its2.ds.richness.summarized$`its2.ds.richnesstabmet$State`,Richness=its2.ds.richness.summarized$`its2.ds.richnesstabmet$Richness`)
library(ggplot2)
library(ggpubr)
#Below I create a boxplot for richness values.
its2.ds.rich<-ggboxplot(its2.ds.richnesstabmet, x="State", y="Richness", outlier.shape=NA)+
geom_jitter(data=its2.ds.richnesstabmet,position=position_jitter(0.2), aes(shape=clone),size=3)+
scale_shape_manual(values=c(16,10,8,7,0,2))+
#geom_text(data=richness.summarized, aes(x=State, y = 25 + Richness, label=Group, fontface="bold"))+
geom_text(data=NULL,aes(x=1, y=220, label="A"))+
geom_text(data=NULL,aes(x=2, y=300, label="B"))+
geom_text(data=NULL,aes(x=3, y=175, label="C"))+
geom_text(data=NULL,aes(x=0.5,y=310,label="Caulosphere: ITS2"),hjust=0)+
theme(axis.line.x=element_blank(),
axis.text.x=element_blank(),
axis.title.x=element_blank(),
axis.ticks.x=element_blank())
its2.ds.rich
Diversity
Observed OTU richness is first calculated for each sample. Then due to the unbalanced nature of our study design, we use a type 3 ANOVA.
#In this chunk we calculate shannon diversity and create a data frame including diversity with metadata
#I first import my rarefied OTU table that contains singletons.
its2.ds.rareotu.wsingles<-read.table("Mothur_output/barkits2otutable.rarefied_jnigra17_fa19_gw.shared", header=TRUE, sep="\t")
#I then remove mothur associated metadata present in the first 3 columns.
its2.ds.rareotu<-as.data.frame(its2.ds.rareotu.wsingles[,4:2697])
#Calculate shannon diversity using the diversity function of vegan
library(vegan)
its2.ds.bac.div<-diversity (its2.ds.rareotu, index="shannon")
#I then create a new data frame to include sample metadata.
its2.ds.diversity.tabmet<-data.frame(State=its2.ds.met.rare$State, clone=its2.ds.met.rare$Clone, sample=its2.ds.met.rare$Group,Shannon=its2.ds.bac.div)
#I first test for a significant interaction between state and clone. Due to the unbalanced design we use a type 3 ANOVA.
library(car)
its2.ds.anova.typ3.div<-aov((Shannon)~State*clone, data=its2.ds.diversity.tabmet)
its2.ds.typ3aov.div<-Anova(its2.ds.anova.typ3.div,type="III")
its2.ds.typ3aov.div## Anova Table (Type III tests)
##
## Response: (Shannon)
## Sum Sq Df F value Pr(>F)
## (Intercept) 21.9749 1 115.9089 5.341e-12 ***
## State 0.1600 2 0.4219 0.6595
## clone 0.3021 4 0.3983 0.8083
## State:clone 0.9347 8 0.6163 0.7574
## Residuals 5.8772 31
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1| Sum Sq | Df | F value | Pr(>F) | |
|---|---|---|---|---|
| (Intercept) | 21.9748957 | 1 | 115.9089223 | 0.0000000 |
| State | 0.1599680 | 2 | 0.4218842 | 0.6595188 |
| clone | 0.3020658 | 4 | 0.3983195 | 0.8082868 |
| State:clone | 0.9347117 | 8 | 0.6162795 | 0.7573833 |
| Residuals | 5.8772159 | 31 | NA | NA |
#The below two lines generate plots to evaluate regression assumptions.
par(mfrow=c(2,2))
plot(its2.ds.anova.typ3.div)
#There was no significant interaction between state and clone and thus we can move on to a type 2 ANOVA.
its2.ds.anova.add.div<-aov((Shannon)~State+clone, data=its2.ds.diversity.tabmet)
its2.ds.add.aov.div<-Anova(its2.ds.anova.add.div,type="III")
its2.ds.add.aov.div## Anova Table (Type III tests)
##
## Response: (Shannon)
## Sum Sq Df F value Pr(>F)
## (Intercept) 41.478 1 237.471 < 2.2e-16 ***
## State 6.802 2 19.473 1.367e-06 ***
## clone 1.087 4 1.556 0.2053
## Residuals 6.812 39
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1| Sum Sq | Df | F value | Pr(>F) | |
|---|---|---|---|---|
| (Intercept) | 41.477903 | 1 | 237.471436 | 0.0000000 |
| State | 6.802322 | 2 | 19.472503 | 0.0000014 |
| clone | 1.087081 | 4 | 1.555953 | 0.2052944 |
| Residuals | 6.811928 | 39 | NA | NA |
#Perform Tukey test
its2.ds.add.aov.div.tukey<-TukeyHSD(its2.ds.anova.add.div)
its2.ds.add.aov.div.tukey.table<-rbind(its2.ds.add.aov.div.tukey$State,its2.ds.add.aov.div.tukey$clone)
its2.ds.add.aov.div.tukey.table## diff lwr upr p adj
## TN-IN 0.696594006 0.3337587 1.0594293 1.004467e-04
## WA-IN -0.162085783 -0.5469308 0.2227592 5.650998e-01
## WA-TN -0.858679788 -1.2215151 -0.4958445 3.237409e-06
## 132-130 0.374186198 -0.2065095 0.9548819 3.647146e-01
## 272-130 0.467950301 -0.1127454 1.0486460 1.652566e-01
## 55-130 0.380377650 -0.2171533 0.9779086 3.768785e-01
## WT-130 0.349082949 -0.1963857 0.8945516 3.715678e-01
## 272-132 0.093764103 -0.4695934 0.6571217 9.891090e-01
## 55-132 0.006191452 -0.5745042 0.5868871 9.999998e-01
## WT-132 -0.025103249 -0.5520760 0.5018695 9.999188e-01
## 55-272 -0.087572651 -0.6682683 0.4931230 9.925080e-01
## WT-272 -0.118867352 -0.6458401 0.4081054 9.665610e-01
## WT-55 -0.031294701 -0.5767633 0.5141739 9.998298e-01| diff | lwr | upr | p adj | |
|---|---|---|---|---|
| TN-IN | 0.6965940 | 0.3337587 | 1.0594293 | 0.0001004 |
| WA-IN | -0.1620858 | -0.5469308 | 0.2227592 | 0.5650998 |
| WA-TN | -0.8586798 | -1.2215151 | -0.4958445 | 0.0000032 |
| 132-130 | 0.3741862 | -0.2065095 | 0.9548819 | 0.3647146 |
| 272-130 | 0.4679503 | -0.1127454 | 1.0486460 | 0.1652566 |
| 55-130 | 0.3803777 | -0.2171533 | 0.9779086 | 0.3768785 |
| WT-130 | 0.3490829 | -0.1963857 | 0.8945516 | 0.3715678 |
| 272-132 | 0.0937641 | -0.4695934 | 0.6571217 | 0.9891090 |
| 55-132 | 0.0061915 | -0.5745042 | 0.5868871 | 0.9999998 |
| WT-132 | -0.0251032 | -0.5520760 | 0.5018695 | 0.9999188 |
| 55-272 | -0.0875727 | -0.6682683 | 0.4931230 | 0.9925080 |
| WT-272 | -0.1188674 | -0.6458401 | 0.4081054 | 0.9665610 |
| WT-55 | -0.0312947 | -0.5767633 | 0.5141739 | 0.9998298 |
#The below two lines generate plots to evaluate regression assumptions.
par(mfrow=c(2,2))
plot(its2.ds.anova.add.div)
#There was no significant interaction between state and clone and thus we can move on to a type 2 ANOVA.
its2.ds.anova.sing.div<-aov((Shannon)~State, data=its2.ds.diversity.tabmet)
its2.ds.sing.aov.div<-Anova(its2.ds.anova.sing.div,type="III")
its2.ds.sing.aov.div## Anova Table (Type III tests)
##
## Response: (Shannon)
## Sum Sq Df F value Pr(>F)
## (Intercept) 120.399 1 655.416 < 2.2e-16 ***
## State 6.810 2 18.535 1.567e-06 ***
## Residuals 7.899 43
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1| Sum Sq | Df | F value | Pr(>F) | |
|---|---|---|---|---|
| (Intercept) | 120.398543 | 1 | 655.41611 | 0.0e+00 |
| State | 6.809521 | 2 | 18.53457 | 1.6e-06 |
| Residuals | 7.899008 | 43 | NA | NA |
#Perform Tukey test
its2.ds.anova.sing.div.tukey<-TukeyHSD(its2.ds.anova.sing.div)
its2.ds.anova.sing.div.tukey.table<-rbind(its2.ds.anova.sing.div.tukey$State)
its2.ds.anova.sing.div.tukey.table## diff lwr upr p adj
## TN-IN 0.6965940 0.3258490 1.0673390 1.225586e-04
## WA-IN -0.1620858 -0.5553202 0.2311486 5.805198e-01
## WA-TN -0.8586798 -1.2294248 -0.4879348 3.823443e-06| diff | lwr | upr | p adj | |
|---|---|---|---|---|
| TN-IN | 0.6965940 | 0.3258490 | 1.0673390 | 0.0001226 |
| WA-IN | -0.1620858 | -0.5553202 | 0.2311486 | 0.5805198 |
| WA-TN | -0.8586798 | -1.2294248 | -0.4879348 | 0.0000038 |
#The below two lines generate plots to evaluate regression assumptions.
par(mfrow=c(2,2))
plot(its2.ds.anova.sing.div)
#I then create a box plot for the richness data by state using ggplot.
library(Hmisc)
library(ggpubr)
its2.ds.diversity.summarized<-summarize(its2.ds.diversity.tabmet$Shannon,by=its2.ds.diversity.tabmet$State, FUN=max)
its2.ds.diversity.summarized<-data.frame(State=its2.ds.diversity.summarized$`its2.ds.diversity.tabmet$State`,Diversity=its2.ds.diversity.summarized$`its2.ds.diversity.tabmet$Shannon`)
library(ggplot2)
library(ggpubr)
#Below I create a boxplot for richness values.
its2.ds.div<-ggboxplot(its2.ds.diversity.tabmet, x="State", y="Shannon", outlier.shape=NA)+
geom_jitter(data=its2.ds.diversity.tabmet,position=position_jitter(0.2), aes(shape=clone),size=3)+
scale_shape_manual(values=c(16,10,8,7,0,2))+
#geom_text(data=diversity.summarized, aes(x=State, y = 0.5 + Diversity, label=group, fontface="bold"))+
geom_text(data=NULL,aes(x=1, y=4, label="A"))+
geom_text(data=NULL,aes(x=2, y=4.4, label="B"))+
geom_text(data=NULL,aes(x=3, y=3.6, label="A"))+
geom_text(data=NULL,aes(x=0.50, y=4.5, label="Caulosphere: ITS2"),hjust=0)+
theme(axis.line.x=element_blank(),
axis.text.x=element_blank(),
axis.title.x=element_blank(),
axis.ticks.x=element_blank())
its2.ds.div
OTU Renaming
To give OTUs more meaningful names, taxonomy strings can be subset to change the OTU ID to unique taxonomic ID. Similar to the relative abundance charts, the taxonomy file must first me added to a transposed OTU table, taxonomy strings need to be split into different classification levels, and the lowest level of classification can be used to replace the original OTU ID. For instance, if an OTU received Ascomycota as its highest level of classification, the OTU ID would be changed to Ascomycota1.
Taxonomy Filtering
Taxonomy is first filtered to include only those OTUs present in the rarefied OTU table with no singletons. Then the OTU table is transposed so that OTU ID is the row name and sample ID is the column name. Taxonomy strings are then added to the OTU table and strings are split by semi-colon.
#Below I am reformatting the OTU table so that I can begin to filter and merge the taxonomy information.
#I first transpose the rarefied OTU table with no singletons.
t.its2.ds.rareotu.nosingles<-as.data.frame(t(its2.ds.rareotu.nosingles))
#The below function is asking specifically for the OTU IDs. These are the row names. I will use these to subset the taxonomy file.
t.its2.ds.rareotu.nosingleslabs<-labels(t.its2.ds.rareotu.nosingles)
its2.ds.taxrare<-subset(its2.taxonomy.ds, rownames(its2.taxonomy.ds) %in% t.its2.ds.rareotu.nosingleslabs[[1]])
its2.ds.taxrareinfo<-(its2.ds.taxrare[,2])
#Below I am using the separate function from tidyr to split the taxonomy strings into columns by semi-colon so that I can rename the OTUs to give them more meaning for downstream analyses.
library(tidyr)
its2.ds.taxonomylabs<-c("kingdom","phylum","class","order","family","genus","species")
t.its2.ds.rareotutabtax<-data.frame(taxonomy=its2.ds.taxrareinfo,t.its2.ds.rareotu.nosingles)
t.its2.ds.rareotutabtaxsep<-separate(t.its2.ds.rareotutabtax,into=its2.ds.taxonomylabs,col=taxonomy,sep=";")Renaming of OTUs
OTUs are then renamed. Prior to official renaming, taxonomy strings are cleaned so that assignment confidence values and taxon leaders (e.g. p__, c__, etc.) are removed.
#Here I am renaming the OTUs using the make.names function. This will make the OTU name something more informative. I am using species to do this.
t.its2.ds.rareotutabtaxsep$species<-sub("\\(.*)","",x=t.its2.ds.rareotutabtaxsep$species)
t.its2.ds.rareotutabtaxsep$species<-sub("(_unclassified)","",x=t.its2.ds.rareotutabtaxsep$species)
t.its2.ds.rareotutabtaxsep$species<-sub("(p__)","",x=t.its2.ds.rareotutabtaxsep$species)
t.its2.ds.rareotutabtaxsep$species<-sub("(c__)","",x=t.its2.ds.rareotutabtaxsep$species)
t.its2.ds.rareotutabtaxsep$species<-sub("(o__)","",x=t.its2.ds.rareotutabtaxsep$species)
t.its2.ds.rareotutabtaxsep$species<-sub("(f__)","",x=t.its2.ds.rareotutabtaxsep$species)
t.its2.ds.rareotutabtaxsep$species<-sub("(g__)","",x=t.its2.ds.rareotutabtaxsep$species)
t.its2.ds.rareotutabtaxsep$species<-sub("(s__)","",x=t.its2.ds.rareotutabtaxsep$species)
t.its2.ds.rareotutabtaxsep$species<-sub("_sp","",x=t.its2.ds.rareotutabtaxsep$species)
t.its2.ds.rareotutabtaxsep$species<-sub("_"," ",x=t.its2.ds.rareotutabtaxsep$species)
rownames(t.its2.ds.rareotutabtaxsep)<-make.names(paste(rownames(t.its2.ds.rareotutabtaxsep),t.its2.ds.rareotutabtaxsep$species,sep="_"))
t.its2.ds.rareotutabtax.rename<-t.its2.ds.rareotutabtaxsep[,8:53]
its2.ds.rareotutabtax.rename<-as.data.frame(t(t.its2.ds.rareotutabtax.rename))library(indicspecies)
set.seed(577)
indicspecies.its2.ds<-multipatt(its2.ds.rareotutabtax.rename, its2.ds.met.rare$State,func="IndVal.g",control=how(nperm=10000))
sink(file="its2.ds.indicators.txt")
summary(indicspecies.its2.ds,indvalcomp=TRUE)##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: IndVal.g
## Significance level (alpha): 0.05
##
## Total number of species: 1578
## Selected number of species: 505
## Number of species associated to 1 group: 464
## Number of species associated to 2 groups: 41
##
## List of species associated to each combination:
##
## Group IN #sps. 119
## A B stat p.value
## Otu0050_Rhinocladiella 0.9968 1.0000 0.998 1e-04 ***
## Otu0070_Ascomycota 0.9958 1.0000 0.998 1e-04 ***
## Otu0023_Diplodia 0.9510 1.0000 0.975 1e-04 ***
## Otu0075_Orbilia 1.0000 0.9286 0.964 1e-04 ***
## Otu0206_Pleosporales 1.0000 0.9286 0.964 1e-04 ***
## Otu0257_Rhinocladiella 0.9349 0.9286 0.932 1e-04 ***
## Otu0071_Pleosporales 1.0000 0.8571 0.926 1e-04 ***
## Otu0086_Ascomycota 0.9767 0.8571 0.915 1e-04 ***
## Otu0080_Ascomycota 0.9661 0.8571 0.910 1e-04 ***
## Otu0121_Candelariaceae 0.9602 0.8571 0.907 1e-04 ***
## Otu0435_Ascomycota 0.9720 0.7857 0.874 1e-04 ***
## Otu0018_Phaeomoniellales 0.8162 0.9286 0.871 1e-04 ***
## Otu0031_Pleosporales 0.8704 0.8571 0.864 1e-04 ***
## Otu0317_Basidiomycota 1.0000 0.7143 0.845 1e-04 ***
## Otu0233_Candelaria.fibrosa 0.9054 0.7857 0.843 1e-04 ***
## Otu0029_Phaeoacremonium 0.9746 0.7143 0.834 0.0002 ***
## Otu0200_Phaeomoniellales 0.8596 0.7857 0.822 1e-04 ***
## Otu0393_Rhinocladiella 0.9300 0.7143 0.815 1e-04 ***
## Otu0169_Physcia.millegrana 0.9149 0.7143 0.808 1e-04 ***
## Otu0125_Lophiostomataceae 0.9895 0.6429 0.798 1e-04 ***
## Otu0210_Pleosporales 0.9829 0.6429 0.795 1e-04 ***
## Otu0280_Ascomycota 0.9651 0.6429 0.788 1e-04 ***
## Otu0342_Ascomycota 0.8282 0.7143 0.769 0.0003 ***
## Otu0160_Phaeomoniellales 0.7447 0.7857 0.765 0.0002 ***
## Otu0084_Herpotrichiellaceae 0.6787 0.8571 0.763 0.0003 ***
## Otu0191_Phaeophyscia 1.0000 0.5714 0.756 0.0002 ***
## Otu0202_Chaetothyriales 1.0000 0.5714 0.756 1e-04 ***
## Otu0587_Nectriaceae 1.0000 0.5714 0.756 1e-04 ***
## Otu0119_Ascomycota 0.9914 0.5714 0.753 1e-04 ***
## Otu0449_Phaeomoniellales 0.8809 0.6429 0.753 1e-04 ***
## Otu0172_Ascomycota 0.9715 0.5714 0.745 0.0004 ***
## Otu0249_Pleosporales 0.9686 0.5714 0.744 1e-04 ***
## Otu0632_Rhinocladiella 0.9340 0.5714 0.731 1e-04 ***
## Otu0189_Phaeomoniellales 0.8065 0.6429 0.720 0.0008 ***
## Otu0461_Herpotrichiellaceae 1.0000 0.5000 0.707 0.0003 ***
## Otu0368_Phaeomoniellales 0.8727 0.5714 0.706 0.0006 ***
## Otu0123_Amphisphaeriaceae 0.9947 0.5000 0.705 0.0002 ***
## Otu0171_Chaetothyriales 0.9923 0.5000 0.704 0.0002 ***
## Otu0505_Rhinocladiella 0.9643 0.5000 0.694 0.0002 ***
## Otu0422_Phaeomoniellales 0.7473 0.6429 0.693 0.0005 ***
## Otu0299_Devriesia.pseudoamericana 0.7347 0.6429 0.687 0.0019 **
## Otu0476_Phaeomoniellales 0.8727 0.5000 0.661 0.0022 **
## Otu0043_Mycena 1.0000 0.4286 0.655 0.0009 ***
## Otu0083_Ascomycota 1.0000 0.4286 0.655 0.0008 ***
## Otu0127_Pleosporales 1.0000 0.4286 0.655 0.0010 ***
## Otu0153_Ascomycota 1.0000 0.4286 0.655 0.0009 ***
## Otu0276_Capnodiales 1.0000 0.4286 0.655 0.0006 ***
## Otu0314_Ascomycota 1.0000 0.4286 0.655 0.0005 ***
## Otu0360_Phaeomoniellales 1.0000 0.4286 0.655 0.0003 ***
## Otu0486_Pleosporales 1.0000 0.4286 0.655 0.0011 **
## Otu0577_Ascomycota 1.0000 0.4286 0.655 0.0010 ***
## Otu0592_Rhinocladiella 1.0000 0.4286 0.655 0.0007 ***
## Otu0593_Tremella 1.0000 0.4286 0.655 0.0007 ***
## Otu0188_Ascomycota 0.9829 0.4286 0.649 0.0027 **
## Otu0272_Pleosporales 0.9682 0.4286 0.644 0.0011 **
## Otu0243_Devriesia 0.9651 0.4286 0.643 0.0016 **
## Otu0298_Sordariomycetes 0.7131 0.5714 0.638 0.0026 **
## Otu0134_Cucurbitariaceae 0.9467 0.4286 0.637 0.0040 **
## Otu0174_Candelariella 0.7994 0.5000 0.632 0.0157 *
## Otu0141_Angustimassarina.acerina 0.9007 0.4286 0.621 0.0205 *
## Otu0090_Ascomycota 0.7373 0.5000 0.607 0.0074 **
## Otu0201_Hypocreales 1.0000 0.3571 0.598 0.0028 **
## Otu0259_Punctelia.rudecta 1.0000 0.3571 0.598 0.0026 **
## Otu0316_Auricularia 1.0000 0.3571 0.598 0.0030 **
## Otu0377_Auriculariales 1.0000 0.3571 0.598 0.0034 **
## Otu0408_Bannozyma 1.0000 0.3571 0.598 0.0037 **
## Otu0585_Ascomycota 1.0000 0.3571 0.598 0.0034 **
## Otu0631_Agaricales 1.0000 0.3571 0.598 0.0031 **
## Otu0638_Capnodiales 1.0000 0.3571 0.598 0.0030 **
## Otu0212_Ascomycota 0.9780 0.3571 0.591 0.0021 **
## Otu0602_Dothideomycetes 0.9536 0.3571 0.584 0.0040 **
## Otu0419_Kockovaella 0.9435 0.3571 0.581 0.0054 **
## Otu0469_Phaeophyscia 0.9435 0.3571 0.581 0.0051 **
## Otu0467_Leptospora 0.9278 0.3571 0.576 0.0069 **
## Otu0333_Pleosporales 0.9224 0.3571 0.574 0.0084 **
## Otu0148_Ceratobasidiaceae 0.8944 0.3571 0.565 0.0455 *
## Otu0207_Herpotrichiellaceae 0.8914 0.3571 0.564 0.0113 *
## Otu0327_Phaeomoniellales 0.7415 0.4286 0.564 0.0120 *
## Otu0197_Fellomyces.mexicanus 0.8196 0.3571 0.541 0.0248 *
## Otu0042_Wickerhamomyces.hampshirensis 1.0000 0.2857 0.535 0.0120 *
## Otu0156_Fusarium.solani 1.0000 0.2857 0.535 0.0124 *
## Otu0328_Ascomycota 1.0000 0.2857 0.535 0.0116 *
## Otu0426_Ceratobasidium 1.0000 0.2857 0.535 0.0130 *
## Otu0454_Orbiliaceae 1.0000 0.2857 0.535 0.0126 *
## Otu0527_Pleosporales 1.0000 0.2857 0.535 0.0135 *
## Otu0555_Lophiostoma.cynaroidis 1.0000 0.2857 0.535 0.0128 *
## Otu0722_Capnodiales 1.0000 0.2857 0.535 0.0126 *
## Otu0740_Ascomycota 1.0000 0.2857 0.535 0.0130 *
## Otu0759_Ascomycota 1.0000 0.2857 0.535 0.0127 *
## Otu0783_Nigrospora.oryzae 1.0000 0.2857 0.535 0.0129 *
## Otu0939_Capnodiales 1.0000 0.2857 0.535 0.0113 *
## Otu0323_Phaeoacremonium 0.9794 0.2857 0.529 0.0084 **
## Otu0999_Phaeosphaeriaceae 0.8571 0.2857 0.495 0.0399 *
## Otu0046_Pleosporales 1.0000 0.2143 0.463 0.0496 *
## Otu0079_Pleosporales 1.0000 0.2143 0.463 0.0445 *
## Otu0457_Auriculariales 1.0000 0.2143 0.463 0.0489 *
## Otu0459_Ceratobasidiaceae 1.0000 0.2143 0.463 0.0499 *
## Otu0472_Ceratobasidium 1.0000 0.2143 0.463 0.0480 *
## Otu0482_Fusicolla.violacea 1.0000 0.2143 0.463 0.0445 *
## Otu0485_Prosthemium.betulinum 1.0000 0.2143 0.463 0.0458 *
## Otu0492_Kockovaella 1.0000 0.2143 0.463 0.0460 *
## Otu0534_Pleosporales 1.0000 0.2143 0.463 0.0480 *
## Otu0590_Pleosporales 1.0000 0.2143 0.463 0.0491 *
## Otu0658_Pleosporales 1.0000 0.2143 0.463 0.0461 *
## Otu0713_Leptospora.rubella 1.0000 0.2143 0.463 0.0483 *
## Otu0778_Ascomycota 1.0000 0.2143 0.463 0.0441 *
## Otu0870_Ceratobasidiaceae 1.0000 0.2143 0.463 0.0457 *
## Otu0936_Pleosporales 1.0000 0.2143 0.463 0.0485 *
## Otu0945_Septobasidiaceae 1.0000 0.2143 0.463 0.0496 *
## Otu0973_Kockovaella 1.0000 0.2143 0.463 0.0490 *
## Otu0980_Ascomycota 1.0000 0.2143 0.463 0.0465 *
## Otu1028_Phaeomoniellales 1.0000 0.2143 0.463 0.0489 *
## Otu1069_Phaeomoniellales 1.0000 0.2143 0.463 0.0480 *
## Otu1099_Myriangium 1.0000 0.2143 0.463 0.0492 *
## Otu1100_Phaeomoniellaceae 1.0000 0.2143 0.463 0.0470 *
## Otu1110_Ascomycota 1.0000 0.2143 0.463 0.0469 *
## Otu1136_Cryptodiscus 1.0000 0.2143 0.463 0.0470 *
## Otu1290_Ascomycota 1.0000 0.2143 0.463 0.0470 *
## Otu1364_Phaeomoniellales 1.0000 0.2143 0.463 0.0481 *
##
## Group TN #sps. 196
## A B stat p.value
## Otu0014_Paraconiothyrium 0.9770 1.0000 0.988 1e-04 ***
## Otu0081_Lecanorales 1.0000 0.9444 0.972 1e-04 ***
## Otu0177_Capnodiales 1.0000 0.9444 0.972 1e-04 ***
## Otu0058_Ascomycota 0.9963 0.9444 0.970 1e-04 ***
## Otu0124_Trichomeriaceae 0.9639 0.9444 0.954 1e-04 ***
## Otu0068_Ascomycota 1.0000 0.8889 0.943 1e-04 ***
## Otu0104_Ascomycota 1.0000 0.8889 0.943 1e-04 ***
## Otu0252_Capnodiales 1.0000 0.8889 0.943 1e-04 ***
## Otu0019_Ascomycota 0.9935 0.8889 0.940 1e-04 ***
## Otu0017_Pleosporales 0.9886 0.8889 0.937 1e-04 ***
## Otu0077_Pleosporales 0.9874 0.8889 0.937 1e-04 ***
## Otu0049_Leotiomycetes 0.9836 0.8889 0.935 1e-04 ***
## Otu0176_Lecanoromycetes 1.0000 0.8333 0.913 1e-04 ***
## Otu0219_Myriangium.citri 1.0000 0.8333 0.913 1e-04 ***
## Otu0015_Ascomycota 0.9995 0.8333 0.913 1e-04 ***
## Otu0167_Ascomycota 0.9946 0.8333 0.910 1e-04 ***
## Otu0163_Endosporium 0.9852 0.8333 0.906 1e-04 ***
## Otu0105_Rachicladosporium 1.0000 0.7778 0.882 1e-04 ***
## Otu0117_Pleosporales 1.0000 0.7778 0.882 1e-04 ***
## Otu0154_Trichomeriaceae 1.0000 0.7778 0.882 1e-04 ***
## Otu0173_Trichomeriaceae 1.0000 0.7778 0.882 1e-04 ***
## Otu0204_Lecanorales 1.0000 0.7778 0.882 1e-04 ***
## Otu0347_Capnodiales 1.0000 0.7778 0.882 1e-04 ***
## Otu0136_Endosporium.aviarium 0.9939 0.7778 0.879 1e-04 ***
## Otu0091_Lecanoromycetes 0.9932 0.7778 0.879 1e-04 ***
## Otu0251_Ascomycota 0.9840 0.7778 0.875 1e-04 ***
## Otu0048_Pleosporales 0.8768 0.8333 0.855 0.0011 **
## Otu0088_Pleosporales 1.0000 0.7222 0.850 1e-04 ***
## Otu0098_Trichomeriaceae 1.0000 0.7222 0.850 1e-04 ***
## Otu0126_Trichomeriaceae 1.0000 0.7222 0.850 1e-04 ***
## Otu0175_Ascomycota 1.0000 0.7222 0.850 1e-04 ***
## Otu0190_Trichomeriaceae 1.0000 0.7222 0.850 1e-04 ***
## Otu0232_Lecanoromycetes 1.0000 0.7222 0.850 1e-04 ***
## Otu0056_Hyperphyscia.adglutinata 0.9954 0.7222 0.848 1e-04 ***
## Otu0270_Strelitziana.africana 0.9868 0.7222 0.844 1e-04 ***
## Otu0195_Capnodiales 0.9694 0.7222 0.837 1e-04 ***
## Otu0215_Lecanoromycetes 1.0000 0.6667 0.816 1e-04 ***
## Otu0236_Lecanoromycetes 1.0000 0.6667 0.816 0.0002 ***
## Otu0303_Chaetothyriales 1.0000 0.6667 0.816 1e-04 ***
## Otu0305_Orbiliales 1.0000 0.6667 0.816 1e-04 ***
## Otu0101_Lophiostoma.fuckelii 0.9934 0.6667 0.814 1e-04 ***
## Otu0269_Ascomycota 0.9901 0.6667 0.812 1e-04 ***
## Otu0238_Capnodiales 0.9901 0.6667 0.812 1e-04 ***
## Otu0423_Neodevriesia.lagerstroemiae 0.9333 0.6667 0.789 1e-04 ***
## Otu0011_Pleosporales 1.0000 0.6111 0.782 1e-04 ***
## Otu0203_Strelitziana.africana 1.0000 0.6111 0.782 1e-04 ***
## Otu0209_Ascomycota 1.0000 0.6111 0.782 1e-04 ***
## Otu0222_Capnodiales 1.0000 0.6111 0.782 1e-04 ***
## Otu0231_Strelitziana.africana 1.0000 0.6111 0.782 1e-04 ***
## Otu0335_Dothideomycetes 1.0000 0.6111 0.782 1e-04 ***
## Otu0355_Basidiomycota 1.0000 0.6111 0.782 1e-04 ***
## Otu0365_Ascomycota 1.0000 0.6111 0.782 1e-04 ***
## Otu0085_Periconia.macrospinosa 0.9128 0.6667 0.780 0.0021 **
## Otu0016_Didymellaceae 0.9937 0.6111 0.779 0.0002 ***
## Otu0100_Herpotrichiellaceae 0.9912 0.6111 0.778 1e-04 ***
## Otu0258_Pleosporales 0.8738 0.6667 0.763 0.0003 ***
## Otu0297_Ascomycota 0.9389 0.6111 0.757 0.0006 ***
## Otu0295_Ascomycota 1.0000 0.5556 0.745 1e-04 ***
## Otu0340_Ascomycota 1.0000 0.5556 0.745 0.0003 ***
## Otu0411_Capnodiales 1.0000 0.5556 0.745 0.0002 ***
## Otu0239_Ascomycota 0.9883 0.5556 0.741 0.0003 ***
## Otu0348_Trichomeriaceae 0.9749 0.5556 0.736 0.0002 ***
## Otu0226_Lecanoromycetes 0.9661 0.5556 0.733 0.0005 ***
## Otu0468_Capnodiales 0.9625 0.5556 0.731 0.0007 ***
## Otu0309_Sphaceloma 0.9561 0.5556 0.729 0.0004 ***
## Otu0234_Sphaceloma 0.8615 0.6111 0.726 0.0009 ***
## Otu0067_Didymellaceae 0.7763 0.6667 0.719 0.0093 **
## Otu0211_Taphrina.inositophila 0.9267 0.5556 0.718 0.0016 **
## Otu0266_Ascomycota 1.0000 0.5000 0.707 0.0003 ***
## Otu0274_Pleosporales 1.0000 0.5000 0.707 0.0002 ***
## Otu0322_Ascomycota 1.0000 0.5000 0.707 0.0006 ***
## Otu0352_Ascomycota 1.0000 0.5000 0.707 0.0003 ***
## Otu0358_Orbiliales 1.0000 0.5000 0.707 1e-04 ***
## Otu0379_Ascomycota 1.0000 0.5000 0.707 0.0003 ***
## Otu0383_Ascomycota 1.0000 0.5000 0.707 0.0002 ***
## Otu0441_Capnodiales 1.0000 0.5000 0.707 1e-04 ***
## Otu0466_Ascomycota 1.0000 0.5000 0.707 1e-04 ***
## Otu0496_Lecanorales 1.0000 0.5000 0.707 0.0005 ***
## Otu0524_Trichomeriaceae 1.0000 0.5000 0.707 0.0003 ***
## Otu0551_Ascomycota 1.0000 0.5000 0.707 0.0004 ***
## Otu0065_Botryosphaeria 0.9973 0.5000 0.706 0.0007 ***
## Otu0116_Diaporthe 0.9931 0.5000 0.705 0.0004 ***
## Otu0273_Ascomycota 0.8814 0.5556 0.700 0.0020 **
## Otu0650_Trichomeriaceae 0.9100 0.5000 0.675 0.0010 ***
## Otu0147_Pleosporales 1.0000 0.4444 0.667 0.0008 ***
## Otu0165_Setomelanomma 1.0000 0.4444 0.667 0.0008 ***
## Otu0291_Lecanorales 1.0000 0.4444 0.667 0.0011 **
## Otu0315_Ascomycota 1.0000 0.4444 0.667 0.0010 ***
## Otu0324_Paraphoma 1.0000 0.4444 0.667 0.0008 ***
## Otu0339_Capnodiales 1.0000 0.4444 0.667 0.0009 ***
## Otu0367_Ascomycota 1.0000 0.4444 0.667 0.0006 ***
## Otu0372_Ascomycota 1.0000 0.4444 0.667 0.0009 ***
## Otu0405_Capnodiales 1.0000 0.4444 0.667 0.0010 ***
## Otu0418_Ascomycota 1.0000 0.4444 0.667 0.0009 ***
## Otu0559_Dothideomycetes 1.0000 0.4444 0.667 0.0012 **
## Otu0582_Myriangium.citri 1.0000 0.4444 0.667 0.0008 ***
## Otu0244_Diaporthe 0.8844 0.5000 0.665 0.0069 **
## Otu0310_Pleosporales 0.8393 0.5000 0.648 0.0040 **
## Otu0607_Ascomycota 0.9333 0.4444 0.644 0.0034 **
## Otu0350_Pleosporales 0.9100 0.4444 0.636 0.0033 **
## Otu0159_Mycoleptodiscus 0.8822 0.4444 0.626 0.0194 *
## Otu0196_Setophoma.chromolaenae 1.0000 0.3889 0.624 0.0028 **
## Otu0312_Ascomycota 1.0000 0.3889 0.624 0.0034 **
## Otu0375_Ascomycota 1.0000 0.3889 0.624 0.0032 **
## Otu0414_Tremella 1.0000 0.3889 0.624 0.0030 **
## Otu0428_Lecanorales 1.0000 0.3889 0.624 0.0020 **
## Otu0451_Pleosporales 1.0000 0.3889 0.624 0.0021 **
## Otu0462_Phaeomoniellales 1.0000 0.3889 0.624 0.0034 **
## Otu0487_Ascomycota 1.0000 0.3889 0.624 0.0016 **
## Otu0576_Ascomycota 1.0000 0.3889 0.624 0.0028 **
## Otu0758_Erythrobasidium 1.0000 0.3889 0.624 0.0020 **
## Otu0767_Ascomycota 1.0000 0.3889 0.624 0.0024 **
## Otu0780_Capnodiales 1.0000 0.3889 0.624 0.0007 ***
## Otu0143_Pestalotiopsis 0.9798 0.3889 0.617 0.0113 *
## Otu0122_Nectriaceae 0.8540 0.4444 0.616 0.0138 *
## Otu0432_Ascomycota 0.9636 0.3889 0.612 0.0058 **
## Otu0187_Phaeosphaeriaceae 0.8974 0.3889 0.591 0.0202 *
## Otu0488_Basidiomycota 0.8861 0.3889 0.587 0.0085 **
## Otu0719_Ascomycota 0.8861 0.3889 0.587 0.0093 **
## Otu0111_Pleosporales 1.0000 0.3333 0.577 0.0124 *
## Otu0225_Pleosporales 1.0000 0.3333 0.577 0.0091 **
## Otu0237_Ascomycota 1.0000 0.3333 0.577 0.0049 **
## Otu0300_Orbilia 1.0000 0.3333 0.577 0.0098 **
## Otu0301_Lecanoromycetes 1.0000 0.3333 0.577 0.0141 *
## Otu0325_Ascomycota 1.0000 0.3333 0.577 0.0107 *
## Otu0386_Trichomeriaceae 1.0000 0.3333 0.577 0.0091 **
## Otu0460_Ascomycota 1.0000 0.3333 0.577 0.0068 **
## Otu0465_Ascomycota 1.0000 0.3333 0.577 0.0068 **
## Otu0501_Lecanorales 1.0000 0.3333 0.577 0.0140 *
## Otu0510_Ascomycota 1.0000 0.3333 0.577 0.0031 **
## Otu0516_Cystobasidiomycetes 1.0000 0.3333 0.577 0.0116 *
## Otu0530_Lecanoromycetes 1.0000 0.3333 0.577 0.0071 **
## Otu0531_Ascomycota 1.0000 0.3333 0.577 0.0076 **
## Otu0595_Ascomycota 1.0000 0.3333 0.577 0.0087 **
## Otu0629_Orbiliales 1.0000 0.3333 0.577 0.0047 **
## Otu0666_Trichomeriaceae 1.0000 0.3333 0.577 0.0052 **
## Otu0668_Pleosporales 1.0000 0.3333 0.577 0.0041 **
## Otu0694_Dothideomycetes 1.0000 0.3333 0.577 0.0086 **
## Otu0698_Trichomeriaceae 1.0000 0.3333 0.577 0.0028 **
## Otu0699_Ascomycota 1.0000 0.3333 0.577 0.0041 **
## Otu0734_Septobasidium 1.0000 0.3333 0.577 0.0070 **
## Otu0942_Trichomeriaceae 1.0000 0.3333 0.577 0.0032 **
## Otu1021_Lecanorales 1.0000 0.3333 0.577 0.0030 **
## Otu0382_Ascomycota 0.9800 0.3333 0.572 0.0128 *
## Otu0025_Cucurbitariaceae 0.9798 0.3333 0.572 0.0117 *
## Otu0436_Trichomeriaceae 0.9561 0.3333 0.565 0.0276 *
## Otu0552_Ascomycota 0.9561 0.3333 0.565 0.0136 *
## Otu0645_Helminthosporium 0.8349 0.3333 0.528 0.0234 *
## Otu0279_Ascomycota 1.0000 0.2778 0.527 0.0236 *
## Otu0374_Chaetothyriales 1.0000 0.2778 0.527 0.0223 *
## Otu0387_Elsinoaceae 1.0000 0.2778 0.527 0.0085 **
## Otu0420_Trichomeriaceae 1.0000 0.2778 0.527 0.0206 *
## Otu0439_Neocucurbitaria 1.0000 0.2778 0.527 0.0199 *
## Otu0448_Sphaceloma 1.0000 0.2778 0.527 0.0104 *
## Otu0515_Ascomycota 1.0000 0.2778 0.527 0.0113 *
## Otu0575_Ascomycota 1.0000 0.2778 0.527 0.0093 **
## Otu0589_Jianyunia.sakaguchii 1.0000 0.2778 0.527 0.0107 *
## Otu0591_Capnodiales 1.0000 0.2778 0.527 0.0092 **
## Otu0609_Dothideomycetes 1.0000 0.2778 0.527 0.0075 **
## Otu0610_Thyridariaceae 1.0000 0.2778 0.527 0.0089 **
## Otu0662_Ascomycota 1.0000 0.2778 0.527 0.0067 **
## Otu0680_Ascomycota 1.0000 0.2778 0.527 0.0090 **
## Otu0688_Chaetothyriales 1.0000 0.2778 0.527 0.0097 **
## Otu0692_Ascomycota 1.0000 0.2778 0.527 0.0092 **
## Otu0702_Ascomycota 1.0000 0.2778 0.527 0.0096 **
## Otu0736_Ascomycota 1.0000 0.2778 0.527 0.0104 *
## Otu0799_Ascomycota 1.0000 0.2778 0.527 0.0088 **
## Otu0854_Lecanorales 1.0000 0.2778 0.527 0.0094 **
## Otu0871_Lecanorales 1.0000 0.2778 0.527 0.0116 *
## Otu0948_Lecanoromycetes 1.0000 0.2778 0.527 0.0087 **
## Otu1048_Pleosporales 1.0000 0.2778 0.527 0.0111 *
## Otu0473_Kurtzmanomyces 0.9511 0.2778 0.514 0.0291 *
## Otu0161_Ascomycota 1.0000 0.2222 0.471 0.0323 *
## Otu0181_Pseudochaetosphaeronema 1.0000 0.2222 0.471 0.0277 *
## Otu0268_Lecanoromycetes 1.0000 0.2222 0.471 0.0334 *
## Otu0453_Spizellomycetales 1.0000 0.2222 0.471 0.0307 *
## Otu0495_Ascomycota 1.0000 0.2222 0.471 0.0346 *
## Otu0506_Lecanoromycetes 1.0000 0.2222 0.471 0.0298 *
## Otu0546_Orbilia.aristata 1.0000 0.2222 0.471 0.0301 *
## Otu0616_Taphrina.inositophila 1.0000 0.2222 0.471 0.0306 *
## Otu0627_Ascomycota 1.0000 0.2222 0.471 0.0344 *
## Otu0721_Orbilia 1.0000 0.2222 0.471 0.0305 *
## Otu0751_Pleosporales 1.0000 0.2222 0.471 0.0280 *
## Otu0761_Trichomeriaceae 1.0000 0.2222 0.471 0.0329 *
## Otu0763_Trichomeriaceae 1.0000 0.2222 0.471 0.0298 *
## Otu0774_Hyperphyscia.adglutinata 1.0000 0.2222 0.471 0.0298 *
## Otu0812_Trichomeriaceae 1.0000 0.2222 0.471 0.0307 *
## Otu0828_Ascomycota 1.0000 0.2222 0.471 0.0317 *
## Otu0841_Ascomycota 1.0000 0.2222 0.471 0.0315 *
## Otu0879_Lecanorales 1.0000 0.2222 0.471 0.0299 *
## Otu0913_Lecanora 1.0000 0.2222 0.471 0.0344 *
## Otu0925_Ascomycota 1.0000 0.2222 0.471 0.0311 *
## Otu0964_Chionosphaeraceae 1.0000 0.2222 0.471 0.0280 *
## Otu0994_Lecanorales 1.0000 0.2222 0.471 0.0296 *
## Otu0995_Ascomycota 1.0000 0.2222 0.471 0.0310 *
## Otu1298_Ascomycota 1.0000 0.2222 0.471 0.0322 *
##
## Group WA #sps. 149
## A B stat p.value
## Otu0007_Pleosporales.fam_Incertae_sedis 1.0000 1.0000 1.000 1e-04 ***
## Otu0008_Melanommataceae 1.0000 1.0000 1.000 1e-04 ***
## Otu0013_Taphrinales 1.0000 1.0000 1.000 1e-04 ***
## Otu0026_Ascomycota 1.0000 1.0000 1.000 1e-04 ***
## Otu0033_Ascomycota 1.0000 1.0000 1.000 1e-04 ***
## Otu0034_Phaeococcomyces 1.0000 1.0000 1.000 1e-04 ***
## Otu0036_Coniothyriaceae 1.0000 1.0000 1.000 1e-04 ***
## Otu0041_Dothideales 1.0000 1.0000 1.000 1e-04 ***
## Otu0057_Aureobasidium.pullulans 1.0000 1.0000 1.000 1e-04 ***
## Otu0059_Filobasidium.wieringae 1.0000 1.0000 1.000 1e-04 ***
## Otu0078_Dothideales 1.0000 1.0000 1.000 1e-04 ***
## Otu0002_Ascomycota 0.9996 1.0000 1.000 1e-04 ***
## Otu0021_Ascomycota 0.9989 1.0000 0.999 1e-04 ***
## Otu0039_Buckleyzyma.aurantiaca 0.9985 1.0000 0.999 1e-04 ***
## Otu0006_Aureobasidium.pullulans 0.9975 1.0000 0.999 1e-04 ***
## Otu0010_Endoconidioma.populi 0.9714 1.0000 0.986 1e-04 ***
## Otu0001_Didymellaceae 0.9647 1.0000 0.982 1e-04 ***
## Otu0044_Ascomycota 1.0000 0.9286 0.964 1e-04 ***
## Otu0095_Filobasidiales 1.0000 0.9286 0.964 1e-04 ***
## Otu0131_Taphrina 1.0000 0.9286 0.964 1e-04 ***
## Otu0326_Taphrinales 1.0000 0.9286 0.964 1e-04 ***
## Otu0009_Alternaria.alternata 0.9126 1.0000 0.955 0.0002 ***
## Otu0194_Sydowia.polyspora 0.9496 0.9286 0.939 1e-04 ***
## Otu0004_Melanommataceae 1.0000 0.8571 0.926 1e-04 ***
## Otu0028_Cryptococcus.cuniculi 1.0000 0.8571 0.926 1e-04 ***
## Otu0054_Microbotryomycetes 1.0000 0.8571 0.926 1e-04 ***
## Otu0074_Kondoa 1.0000 0.8571 0.926 1e-04 ***
## Otu0110_Vishniacozyma.dimennae 1.0000 0.8571 0.926 1e-04 ***
## Otu0118_Knufia 1.0000 0.8571 0.926 1e-04 ***
## Otu0146_Gelidatrema 1.0000 0.8571 0.926 1e-04 ***
## Otu0157_Dothideales 1.0000 0.8571 0.926 1e-04 ***
## Otu0198_Phaeococcomyces 1.0000 0.8571 0.926 1e-04 ***
## Otu0329_Ascomycota 1.0000 0.8571 0.926 1e-04 ***
## Otu0032_Phaeosphaeriaceae 0.9980 0.8571 0.925 1e-04 ***
## Otu0112_Filobasidium.magnum 0.9424 0.8571 0.899 1e-04 ***
## Otu0094_Alternaria.metachromatica 1.0000 0.7857 0.886 1e-04 ***
## Otu0108_Orbilia 1.0000 0.7857 0.886 1e-04 ***
## Otu0138_Dothideomycetes 1.0000 0.7857 0.886 1e-04 ***
## Otu0162_Cystobasidiomycetes 1.0000 0.7857 0.886 1e-04 ***
## Otu0205_Ascomycota 1.0000 0.7857 0.886 1e-04 ***
## Otu0332_Genolevuria 1.0000 0.7857 0.886 1e-04 ***
## Otu0076_Endoconidioma.populi 0.9944 0.7857 0.884 1e-04 ***
## Otu0264_Didymellaceae 0.9118 0.7857 0.846 1e-04 ***
## Otu0120_Leucosporidiales 1.0000 0.7143 0.845 1e-04 ***
## Otu0140_Camarosporidiella 1.0000 0.7143 0.845 1e-04 ***
## Otu0170_Cryptococcus.cuniculi 1.0000 0.7143 0.845 1e-04 ***
## Otu0182_Taphrina.carpini 1.0000 0.7143 0.845 1e-04 ***
## Otu0289_Ascomycota 0.9643 0.7143 0.830 1e-04 ***
## Otu0035_Microbotryomycetes 1.0000 0.6429 0.802 1e-04 ***
## Otu0097_Phaeosphaeriaceae 1.0000 0.6429 0.802 1e-04 ***
## Otu0149_Phaeomoniellales 1.0000 0.6429 0.802 1e-04 ***
## Otu0185_Cystobasidiomycetes 1.0000 0.6429 0.802 1e-04 ***
## Otu0250_Chionosphaeraceae 1.0000 0.6429 0.802 1e-04 ***
## Otu0282_Ascomycota 1.0000 0.6429 0.802 0.0002 ***
## Otu0286_Ramimonilia.apicalis 1.0000 0.6429 0.802 1e-04 ***
## Otu0137_Alternaria.subcucurbitae 1.0000 0.5714 0.756 1e-04 ***
## Otu0213_Didymellaceae 1.0000 0.5714 0.756 1e-04 ***
## Otu0235_Basidiomycota 1.0000 0.5714 0.756 1e-04 ***
## Otu0240_Vishniacozyma.victoriae 1.0000 0.5714 0.756 1e-04 ***
## Otu0313_Ascomycota 1.0000 0.5714 0.756 0.0002 ***
## Otu0349_Ascomycota 1.0000 0.5714 0.756 0.0002 ***
## Otu0024_Chaetothyriales 0.8969 0.5714 0.716 0.0100 **
## Otu0047_Wickerhamomyces.hampshirensis 1.0000 0.5000 0.707 1e-04 ***
## Otu0073_Ustilago.hordei 1.0000 0.5000 0.707 0.0003 ***
## Otu0115_Geosmithia 1.0000 0.5000 0.707 1e-04 ***
## Otu0168_Erythrobasidiales 1.0000 0.5000 0.707 1e-04 ***
## Otu0217_Vishniacozyma.dimennae 1.0000 0.5000 0.707 0.0002 ***
## Otu0246_Tremellales 1.0000 0.5000 0.707 0.0002 ***
## Otu0344_Filobasidium.globisporum 1.0000 0.5000 0.707 1e-04 ***
## Otu0399_Ascomycota 1.0000 0.5000 0.707 0.0002 ***
## Otu0517_Cryptococcus.cuniculi 1.0000 0.5000 0.707 0.0002 ***
## Otu0106_Amphisphaeriaceae 0.9616 0.5000 0.693 0.0083 **
## Otu0051_Pleosporales 1.0000 0.4286 0.655 0.0007 ***
## Otu0184_Corticifraga.peltigerae 1.0000 0.4286 0.655 0.0011 **
## Otu0186_Comoclathris.rosae 1.0000 0.4286 0.655 0.0012 **
## Otu0192_Leptosphaeriaceae 1.0000 0.4286 0.655 0.0005 ***
## Otu0229_Stemphylium 1.0000 0.4286 0.655 0.0012 **
## Otu0256_Microbotryomycetes 1.0000 0.4286 0.655 0.0005 ***
## Otu0346_Pleosporales 1.0000 0.4286 0.655 0.0010 ***
## Otu0362_Ascomycota 1.0000 0.4286 0.655 0.0005 ***
## Otu0364_Lecanorales 1.0000 0.4286 0.655 0.0010 ***
## Otu0376_Dioszegia 1.0000 0.4286 0.655 0.0004 ***
## Otu0378_Udeniomyces.puniceus 1.0000 0.4286 0.655 0.0007 ***
## Otu0427_Taphrinales 1.0000 0.4286 0.655 0.0010 ***
## Otu0474_Pleosporales 1.0000 0.4286 0.655 0.0006 ***
## Otu0513_Tremellales 1.0000 0.4286 0.655 0.0007 ***
## Otu0514_Basidiomycota 1.0000 0.4286 0.655 0.0010 ***
## Otu0543_Teloschistaceae 1.0000 0.4286 0.655 0.0005 ***
## Otu0053_Leotiomycetes 1.0000 0.3571 0.598 0.0039 **
## Otu0055_Melanommataceae 1.0000 0.3571 0.598 0.0028 **
## Otu0099_Leptosphaeriaceae 1.0000 0.3571 0.598 0.0041 **
## Otu0152_Pleosporales 1.0000 0.3571 0.598 0.0027 **
## Otu0224_Dematiopleospora 1.0000 0.3571 0.598 0.0031 **
## Otu0278_Basidiomycota 1.0000 0.3571 0.598 0.0024 **
## Otu0308_Cystobasidiomycetes 1.0000 0.3571 0.598 0.0039 **
## Otu0319_Ascomycota 1.0000 0.3571 0.598 0.0029 **
## Otu0331_Knufia 1.0000 0.3571 0.598 0.0037 **
## Otu0359_Orbiliaceae 1.0000 0.3571 0.598 0.0031 **
## Otu0370_Dothideales 1.0000 0.3571 0.598 0.0032 **
## Otu0381_Vishniacozyma.dimennae 1.0000 0.3571 0.598 0.0033 **
## Otu0401_Cyrenella.elegans 1.0000 0.3571 0.598 0.0029 **
## Otu0412_Naganishia.albida 1.0000 0.3571 0.598 0.0026 **
## Otu0463_Cystobasidiomycetes 1.0000 0.3571 0.598 0.0033 **
## Otu0497_Dioszegia 1.0000 0.3571 0.598 0.0033 **
## Otu0536_Herpotrichia 1.0000 0.3571 0.598 0.0025 **
## Otu0545_Taphrinales 1.0000 0.3571 0.598 0.0035 **
## Otu0558_Knufia 1.0000 0.3571 0.598 0.0032 **
## Otu0565_Basidiomycota 1.0000 0.3571 0.598 0.0031 **
## Otu0151_Chaetosphaeronema 1.0000 0.2857 0.535 0.0129 *
## Otu0179_Phaeomoniellaceae 1.0000 0.2857 0.535 0.0124 *
## Otu0248_Ascomycota 1.0000 0.2857 0.535 0.0113 *
## Otu0304_Leptosphaeria.rubefaciens 1.0000 0.2857 0.535 0.0134 *
## Otu0400_Cystobasidiomycetes 1.0000 0.2857 0.535 0.0128 *
## Otu0440_Knufia 1.0000 0.2857 0.535 0.0110 *
## Otu0444_Pleosporales 1.0000 0.2857 0.535 0.0109 *
## Otu0509_Melanommataceae 1.0000 0.2857 0.535 0.0121 *
## Otu0511_Tremellales 1.0000 0.2857 0.535 0.0114 *
## Otu0522_Didymellaceae 1.0000 0.2857 0.535 0.0120 *
## Otu0539_Dioszegia 1.0000 0.2857 0.535 0.0129 *
## Otu0540_Melanommataceae 1.0000 0.2857 0.535 0.0145 *
## Otu0730_Ascomycota 1.0000 0.2857 0.535 0.0112 *
## Otu0735_Ascomycota 1.0000 0.2857 0.535 0.0130 *
## Otu0769_Melanommataceae 1.0000 0.2857 0.535 0.0124 *
## Otu0796_Didymellaceae 1.0000 0.2857 0.535 0.0115 *
## Otu0843_Taphrinales 1.0000 0.2857 0.535 0.0124 *
## Otu0281_Leptosphaeriaceae 1.0000 0.2143 0.463 0.0452 *
## Otu0288_Thyridariaceae 1.0000 0.2143 0.463 0.0479 *
## Otu0290_Cadophora 1.0000 0.2143 0.463 0.0449 *
## Otu0398_Pleurophoma 1.0000 0.2143 0.463 0.0481 *
## Otu0417_Cystofilobasidium.capitatum 1.0000 0.2143 0.463 0.0479 *
## Otu0424_Microbotryomycetes 1.0000 0.2143 0.463 0.0468 *
## Otu0425_Ascomycota 1.0000 0.2143 0.463 0.0484 *
## Otu0573_Microbotryomycetes 1.0000 0.2143 0.463 0.0481 *
## Otu0578_Microbotryomycetes 1.0000 0.2143 0.463 0.0499 *
## Otu0598_Microbotryomycetes 1.0000 0.2143 0.463 0.0486 *
## Otu0611_Melanommataceae 1.0000 0.2143 0.463 0.0469 *
## Otu0612_Ascomycota 1.0000 0.2143 0.463 0.0498 *
## Otu0613_Chaetothyriales 1.0000 0.2143 0.463 0.0467 *
## Otu0682_Cystobasidiomycetes 1.0000 0.2143 0.463 0.0478 *
## Otu0711_Aureobasidium.pullulans 1.0000 0.2143 0.463 0.0498 *
## Otu0714_Microbotryomycetes 1.0000 0.2143 0.463 0.0448 *
## Otu0729_Pleosporales.fam_Incertae_sedis 1.0000 0.2143 0.463 0.0472 *
## Otu0744_Melanommataceae 1.0000 0.2143 0.463 0.0484 *
## Otu0786_Basidiomycota 1.0000 0.2143 0.463 0.0464 *
## Otu0787_Ascomycota 1.0000 0.2143 0.463 0.0495 *
## Otu0809_Cystobasidiomycetes 1.0000 0.2143 0.463 0.0453 *
## Otu1090_Mycosphaerella.tassiana 1.0000 0.2143 0.463 0.0473 *
## Otu1140_Pleosporales 1.0000 0.2143 0.463 0.0474 *
## Otu2224_Endoconidioma.populi 1.0000 0.2143 0.463 0.0488 *
##
## Group IN+TN #sps. 41
## A B stat p.value
## Otu0003_Phaeomoniellales 1.0000 1.0000 1.000 1e-04 ***
## Otu0022_Helminthosporium.asterinum 1.0000 1.0000 1.000 1e-04 ***
## Otu0027_Trichomeriaceae 1.0000 1.0000 1.000 1e-04 ***
## Otu0040_Helminthosporium 1.0000 0.9688 0.984 1e-04 ***
## Otu0005_Phaeomoniellales 1.0000 0.9375 0.968 1e-04 ***
## Otu0069_Didymosphaeriaceae 1.0000 0.9375 0.968 1e-04 ***
## Otu0030_Trichomeriaceae 0.9991 0.9375 0.968 1e-04 ***
## Otu0103_Trichomeriaceae 1.0000 0.8438 0.919 1e-04 ***
## Otu0064_Arthrocatena.tenebrio 0.9966 0.8438 0.917 0.0002 ***
## Otu0063_Pleosporales 1.0000 0.7812 0.884 1e-04 ***
## Otu0082_Physcia 1.0000 0.7812 0.884 0.0002 ***
## Otu0135_Ascomycota 1.0000 0.7500 0.866 1e-04 ***
## Otu0061_Ascomycota 0.9980 0.7500 0.865 1e-04 ***
## Otu0066_Rhinocladiella 0.9979 0.7500 0.865 1e-04 ***
## Otu0130_Rhinocladiella 1.0000 0.7188 0.848 1e-04 ***
## Otu0128_Ascomycota 1.0000 0.6875 0.829 1e-04 ***
## Otu0230_Capnodiales 1.0000 0.6562 0.810 0.0002 ***
## Otu0093_Phaeophyscia 0.9792 0.6562 0.802 0.0011 **
## Otu0113_Teichosporaceae 1.0000 0.6250 0.791 0.0012 **
## Otu0183_Phaeomoniellales 1.0000 0.5938 0.771 0.0003 ***
## Otu0038_Orbilia 1.0000 0.5625 0.750 0.0032 **
## Otu0216_Didymellaceae 0.9809 0.5625 0.743 0.0036 **
## Otu0114_Xylariales 1.0000 0.5312 0.729 0.0034 **
## Otu0129_Microcera.rubra 1.0000 0.5312 0.729 0.0065 **
## Otu0221_Physciella.chloantha 1.0000 0.5312 0.729 0.0028 **
## Otu0096_Pleosporales 1.0000 0.5000 0.707 0.0049 **
## Otu0155_Tubeufia.cerea 1.0000 0.5000 0.707 0.0085 **
## Otu0037_Caliciopsis.valentina 1.0000 0.4688 0.685 0.0105 *
## Otu0062_Pleosporales 1.0000 0.4688 0.685 0.0116 *
## Otu0265_Dothideomycetes 1.0000 0.4688 0.685 0.0082 **
## Otu0369_Ascomycota 1.0000 0.4688 0.685 0.0049 **
## Otu0532_Capnodiales 1.0000 0.4688 0.685 0.0056 **
## Otu0166_Paraphoma 1.0000 0.4375 0.661 0.0215 *
## Otu0357_Phaeomoniellaceae 1.0000 0.4375 0.661 0.0086 **
## Otu0330_Ascomycota 1.0000 0.4062 0.637 0.0140 *
## Otu0343_Helminthosporium.asterinum 1.0000 0.4062 0.637 0.0170 *
## Otu0445_Bispora.betulina 1.0000 0.3750 0.612 0.0253 *
## Otu0223_Ascomycota 1.0000 0.3438 0.586 0.0499 *
## Otu0431_Orbiliales 1.0000 0.3438 0.586 0.0358 *
## Otu0458_Chionosphaeraceae 1.0000 0.3438 0.586 0.0380 *
## Otu0563_Chionosphaeraceae 1.0000 0.3438 0.586 0.0417 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## State Clone Group
## 1 IN 272 IN_MCB10_272s_S46_L001_R1_001
## 2 IN 272 IN_MCB11_272s_S7_L001_R1_001
## 3 IN 132 IN_MCB16_132s_S70_L001_R1_001
## 4 IN 132 IN_MCB17_132s_S78_L001_R1_001
## 5 IN 272 IN_MCB24_272s_S94_L001_R1_001
## 6 IN 130 IN_MCB26_130s_S62_L001_R1_001
## 7 IN 132 IN_MCB27_132s_S86_L001_R1_001
## 8 IN WT IN_MCB28_WTs_S15_L001_R1_001
## 9 IN 130 IN_MCB2_130s_S38_L001_R1_001
## 10 IN 55 IN_MCB6_55s_S14_L001_R1_001
## 11 IN 55 IN_MCB7_55s_S22_L001_R1_001
## 12 IN 130 IN_MCB9_130s_S54_L001_R1_001
## 13 IN WT IN_WT_MCB_29s_S23_L001_R1
## 14 IN WT IN_WT_MCB_33s_S31_L001_R1
## 15 TN 130 TN_130_As_S71_L001_R1_001
## 16 TN 130 TN_130_Bs_S79_L001_R1_001
## 17 TN 130 TN_130_Cs_S87_L001_R1_001
## 18 TN 132 TN_132_As_S95_L001_R1_001
## 19 TN 132 TN_132_Bs_S47_L001_R1_001
## 20 TN 132 TN_132_Cs_S8_L001_R1_001
## 21 TN 272 TN_272_As_S16_L001_R1_001
## 22 TN 272 TN_272_Bs_S24_L001_R1_001
## 23 TN 272 TN_272_Cs_S32_L001_R1_001
## 24 TN 55 TN_55_As_S39_L001_R1_001
## 25 TN 55 TN_55_Bs_S55_L001_R1_001
## 26 TN 55 TN_55_Cs_S63_L001_R1_001
## 27 TN WT TN_LS_1s_S40_L001_R1_001
## 28 TN WT TN_LS_2s_S56_L001_R1_001
## 29 TN WT TN_LS_3s_S64_L001_R1_001
## 30 TN WT TN_WT_MB_19s_S72_L001_R1
## 31 TN WT TN_WT_MB_20s_S80_L001_R1
## 32 TN WT TN_WT_MB_21s_S88_L001_R1
## 33 WA 272 WA_BNL17_272s_S77_L001_R1_001
## 34 WA 272 WA_BNL18_272s_S85_L001_R1_001
## 35 WA 55 WA_BNL19_55s_S13_L001_R1_001
## 36 WA 130 WA_BNL20_130s_S29_L001_R1_001
## 37 WA WT WA_BNL21_WTs_S93_L001_R1_001
## 38 WA WT WA_BNL22_WTs_S45_L001_R1_001
## 39 WA WT WA_BNL23_WTs_S6_L001_R1_001
## 40 WA 272 WA_RN10_272s_S69_L001_R1_001
## 41 WA 55 WA_RN1_55s_S44_L001_R1_001
## 42 WA 55 WA_RN2_55s_S5_L001_R1_001
## 43 WA 130 WA_RN4_130s_S21_L001_R1_001
## 44 WA 132 WA_RN7_132s_S37_L001_R1_001
## 45 WA 132 WA_RN8_132s_S53_L001_R1_001
## 46 WA 132 WA_RN9_132s_S61_L001_R1_001tcd.status<-c("Negative","Negative","Negative","Negative","Negative","Positive","Positive","Positive","Positive","Negative","Positive","Positive","Positive","Positive")
set.seed(577)
indicspecies.its2.disease<-multipatt(its2.ds.rareotutabtax.rename[33:46,], tcd.status,func="IndVal.g",control=how(nperm=10000))
sink(file="its2.ds.tcdindicators.txt")
summary(indicspecies.its2.disease,indvalcomp=TRUE)##
## Multilevel pattern analysis
## ---------------------------
##
## Association function: IndVal.g
## Significance level (alpha): 0.05
##
## Total number of species: 1578
## Selected number of species: 16
## Number of species associated to 1 group: 16
##
## List of species associated to each combination:
##
## Group Negative #sps. 7
## A B stat p.value
## Otu0349_Ascomycota 0.9256 1.0000 0.962 0.0023 **
## Otu0184_Corticifraga.peltigerae 0.9941 0.8333 0.910 0.0073 **
## Otu0108_Orbilia 0.8101 1.0000 0.900 0.0273 *
## Otu0182_Taphrina.carpini 0.7864 1.0000 0.887 0.0334 *
## Otu0151_Chaetosphaeronema 1.0000 0.6667 0.816 0.0149 *
## Otu0359_Orbiliaceae 0.9730 0.6667 0.805 0.0258 *
## Otu0319_Ascomycota 0.9552 0.6667 0.798 0.0158 *
##
## Group Positive #sps. 9
## A B stat p.value
## Otu0035_Microbotryomycetes 0.9720 1.0000 0.986 0.0013 **
## Otu0185_Cystobasidiomycetes 0.9692 1.0000 0.984 0.0014 **
## Otu0074_Kondoa 0.9261 1.0000 0.962 0.0022 **
## Otu0047_Wickerhamomyces.hampshirensis 1.0000 0.8750 0.935 0.0043 **
## Otu0115_Geosmithia 1.0000 0.8750 0.935 0.0043 **
## Otu0094_Alternaria.metachromatica 0.8467 1.0000 0.920 0.0113 *
## Otu0192_Leptosphaeriaceae 1.0000 0.7500 0.866 0.0183 *
## Otu0256_Microbotryomycetes 1.0000 0.7500 0.866 0.0206 *
## Otu0474_Pleosporales 1.0000 0.7500 0.866 0.0102 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1its2.ds.indicators<-read.table("IndicatorAnalysis/caulosphere_indicatoranalysis_sp20_ajo.csv",header=TRUE,sep=",")
library(tidyr)
its2.ds.indicators.long<-as.data.frame(pivot_longer(its2.ds.indicators, cols=c("Specificity","Sensitivity","Indicator.Value"),"Statistic"))
its2.ds.indicators.long$OTU<-sub("\\."," ",its2.ds.indicators.long$OTU)
in.its2.ds.indicators<-subset(its2.ds.indicators.long,its2.ds.indicators.long[,3]=="IN")
tn.its2.ds.indicators<-subset(its2.ds.indicators.long,its2.ds.indicators.long[,3]=="TN")
tn.in.its2.ds.indicators<-subset(its2.ds.indicators.long,its2.ds.indicators.long[,3]=="IN+TN")
wa.its2.ds.indicators<-subset(its2.ds.indicators.long,its2.ds.indicators.long[,3]=="WA")
tcd.positive.its2.ds.indicators<-subset(its2.ds.indicators.long,its2.ds.indicators.long[,3]=="Positive")
tcd.negative.its2.ds.indicators<-subset(its2.ds.indicators.long,its2.ds.indicators.long[,3]=="Negative")
its2.ds.indicator.mat<-its2.ds.indicators[c(1:10,21:30,41:50,61:70,81:94),2:4]
rownames(its2.ds.indicator.mat)<-its2.ds.indicators$OTU[c(1:10,21:30,41:50,61:70,81:94)]
colnames(its2.ds.indicator.mat)<-colnames(its2.ds.indicators)[2:4]
its2.ds.indicator.row.dat<-data.frame(State=its2.ds.indicators$Group,Mycoparasite=as.factor(its2.ds.indicators$Mycoparasite),Plant.Pathogen=as.factor(its2.ds.indicators$Plant.Pathogen), Saprotroph=as.factor(its2.ds.indicators$Saprotroph))
rownames(its2.ds.indicator.row.dat)<-make.names(its2.ds.indicators$OTU)
annotation_colors<-list(
Saprotroph=c(Yes="dodgerblue4", No="white"),
Plant.Pathogen = c(Yes="maroon", No="white"),
Mycoparasite = c(Yes="khaki4", No="white"),
State=c(IN="springgreen3", IN.TN="slateblue2", TN="skyblue1", WA="violet", Positive="darkslategrey",Negative="olivedrab1")
)
my_colors<-colorRampPalette(colors=c("darkblue","lightblue"))
library(ClassDiscovery)## Loading required package: cluster## Loading required package: oompaBaselibrary(pheatmap)
its.caulosphere.heatmap<-pheatmap(its2.ds.indicator.mat,
annotation_row=its2.ds.indicator.row.dat,
color=gray.colors(15),
cluster_rows=FALSE,
cluster_cols=FALSE,
gaps_row=c(10,20,30,40,49),
annotation_colors=annotation_colors,
cellwidth=10)